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Inversion  of  Waveforms  for  Extreme  Source  Models  with  an 
Application  to  the  Isotropic  Moment  Tensor  Component 

D.  W.  Vasco  * 

L.  R.  Johnson 

Center  for  Computational  Seismology 
Lawrence  Berkeley  Laboratory 
and 

Department  of  Geology  and  Geophysics 
University  of  California 
Berkeley,  CA  94720 

Abstract 

We  have  developed  a  new  approach  to  the  inversion  of  waveform  data  for  the  time- varying 
moment  tensor.  The  method  produces  the  source  model  which  minimizes  the  modulus 
squared  of  any  linear  combination  of  moment  tensor  components,  subject  to  the  constraint 
that  the  data  are  satisfied  within  specified  confidence  intervals.  This  method  allows  the 
determination  of  possible  source  models  other  than  the  least  squares  solution,  enabling  one 
to  determine  the  significance  of  certain  moment  tensor  properties,  for  example,  the  presence 
or  absence  of  a  volume  change  (isotropic  component)  in  the  source.  Synthetic  te'Ts  were 
used  to  examine  the  effect  of  microseismic  noise  and  lateral  heterogeneity  on  the  extreme 
models  of  the  isotropic  component.  Lateral  heterogeneity  is  found  o  have  a  strong  effect 
on  the  estimation  of  the  isotropic  component  of  the  moment  tensor. 

The  method  was  tested  by  using  long-period  waveforms  from  the  Global  Digital  Seis¬ 
mic  Network  to  estimate  the  isotropic  part  of  the  moment  tensor  of  a  deep  Bonin  Islands 
earthquake.  Modelling  indicates  that  more  than  10%  of  the  mechanism  might  have  to  be 
isotropic  for  detection  of  volume  change  in  the  presence  of  10%  random  noise  and  only  2% 
lateral  heterogeneity.  The  least-squares  solution  indicates  that  a  relatively  large  change  in 
volume  was  involved  in  the  source  mechanism.  However,  the  minimum  extreme  solution 
shows  that  this  volume  change  is  not  actually  required  by  the  data  and  thus  may  not  be 
significant.  The  method  was  also  tested  on  near  source  data  from  the  nuclear  explosion 
Harzer.  In  this  case,  in  spite  of  fairly  large  error  bounds,  it  can  be  concluded  that  the 
source  has  a  clear  explosive  component.  ^ 

*Now  at: 

Earth  Sciences  Division 

Air  Force  Geophysics  Laboratory 

Hanscom  Air  Force  Base,  MA  01731 
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INTRODUCTION 


A  major  problem  in  seismology  is  the  determination  of  the  nature  of  seismic  sources.  The 
analysis  of  earthquake  waveform  data  is  the  chief  method  of  studying  the  faulting  process. 
Other  methods  such  as  deep  drilling  and  surface  strain  measurements  are  more  expensive 
and  time  consuming.  Furthermore,  they  do  not  offer  a  dynamic  picture  of  the  faulting 
process.  Therefore  one  must  turn  to  the  elastic  wave  radiation  for  details  of  the  faulting 
process.  One  portion  of  the  seismic  spectrum  in  particular,  that  occupied  by  body  waves, 
will  be  examined  in  this  study  of  seismic  sources.  We  have  chosen  body  waveforms  because  of 
the  higher  resolution  they  contain  and  because  of  the  high  quality  digital  waveforms  which 
are  becoming  available  from  global  digital  networks  such  as  the  Global  Digital  Seismic 
Network  (GDSN).  Waveform  inversion  makes  use  of  the  whole  seismograms,  using  all  the 
information  contained  in  the  time  series. 

With  few  exceptions,  most  moment  tensor  inversions  of  waveforms  have  relied  on  some 
form  of  least  squares  to  derive  a  "best  fitting”  solution  (Gilbert  and  Dziewonski  1975,  Gilbert 
and  Buland  1976,  McCowen  1976,  Kanamori  and  Given  1981,  Stump  and  Johnson  1977, 
Dziewonski  et  al.  1981,  and  Sipkin  1982).  The  notable  exceptions  to  this  approach  are  the 
Lx  minimization  of  Fitch  et  al.  (1980)  and  Tanimoto  and  Kanamori  (1986)  and  the  linear 
programming  approach  of  Julian  (1986).  In  an  approach  related  to  the  latter  paper  we  put 
forth  an  inversion  method  to  examine  the  range  of  time- varying  moment  tensor  models  which 
agree,  to  within  specified  confidence  bounds,  with  observed  body-waveform  data.  Using  this 
technique  it  is  possible  to  derive  extreme  models,  that  is,  models  which  make  some  property 
of  the  source  a  minimum  or  maximum.  As  will  be  detailed  later,  the  properties  will  be  in 
terms  of  the  modulus  squared  of  linear  combinations  of  the  moment  tensor  components. 
Many  properties  can  be  put  in  the  form  of  a  linear  combination  of  moment  components. 
Six  examples  are  given  in  Julian  (1986),  these  include  the  explosiveness,  thrust-like  nature, 
horizontal  extensiveness,  and  vertical  compensated  linear  vector  dipole  nature  of  the  source. 
These  properties  can  be  used  to  answer  interesting  geophysical  questions,  for  example,  to 
decide  if  an  event  is  a  nuclear  explosion  rather  than  an  earthquake. 

One  new  feature  of  our  approach  is  the  use  of  quadratic  programming  for  minimizing  the 
modulus  squared  of  the  linear  combinations  of  moment  tensor  components  rather  than  us¬ 
ing  linear  programming  for  minimizing  a  linear  combination  of  the  components.  Quadratic 
programming  is  the  minimization  or  maximization  of  a  quadratic  functional  subject  to  lin¬ 
ear  equality  and  inequality  constraints  (Dantzig  1963).  This  approach  was  taken  because 
of  problems  in  extending  the  linear  programming  approach  to  waveform  inversion.  Futher- 
more,  it  is  difficult  to  interpret  linear  combinations  of  moment  tensor  components  in  the 
complex  frequency  domain  and  to  transform  these  properties  into  the  time  domain.  As  will 
be  shown,  through  the  use  of  Parseval’s  theorem,  properties  derved  through  the  quadratic 
programming  approach  may  be  conveniently  interpreted  in  both  the  frequency  domain  and 
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the  time  domain.  In  addition,  the  quadratic  programming  method  takes  the  same  order 
of  time  as  the  generalized  least  squares  approach  and  has  proven  to  be  very  efficient  for 
waveform  inversion.  We  explicitly  solve  for  the  time- varying  moment  tensor,  not  assuming 
a  source  time  function.  The  use  of  the  time-varying  moment  tensor  is  a  key  feature  of  the 
method  because  it  allows  for  motion  over  curved  and  complicated  fault  surfaces.  Multiple 
rupture  events  also  pose  no  difficulty  in  the  time- varying  formulation. 

An  application  of  extreme  models,  the  one  principally  addressed  in  this  paper,  is  to 
determine  the  significance  of  the  isotropic  component  in  the  moment  tensor.  This  has 
importance  in  at  least  two  areas:  nuclear  explosion  seismology  and  the  determination  of 
deep  and  intermediate  earthquake  mechanisms.  Nuclear  explosions  are  known  to  be  mainly 
dilatational  sources.  However,  because  of  source  complexity  and  scattering  due  to  lateral 
heterogeneities,  the  moment  tensor  contains  non-dilatational  components.  It  is  necessary  to 
determine  the  significance  of  these  other  components.  Furthermore,  extreme  models  allow 
the  determination  of  the  smallest  and  the  largest  explosive  solutions.  This  is  useful  for  the 
estimation  of  bounds  on  yields. 

In  general,  earthquakes  are  thought  to  have  double  couple  mechanisms.  Occasionally 
however,  investigators  describe  events  which  contain  non-double  couple  components  such 
as  an  isotropic  trace  (Dziewonski  and  Gilbert  1974,  Silver  and  Jordan  1982).  In  particular, 
deep  and  intermediate  earthquakes  in  subduction  zones  may  be  associated  with  volume 
decrease.  The  very  high  pressures  and  the  possibility  of  phase  changes  make  the  fault 
mechanics  very  difficult  to  model.  Therefore,  the  nature  of  the  faulting  deep  within  the  earth 
is  still  an  open  question.  The  role  of  volume  change  can  be  addressed  by  the  calculation 
of  extreme  models.  For  example,  it  is  possible  to  compute  the  model  with  the  minimum 
total  squared  trace  amplitude  and  to  compare  the  trace  of  this  model  with  the  deviatoric 
component.  This  gives  an  indication  of  the  relative  volume  change  associated  with  the 
event.  The  computation  of  such  a  model  results  in  a  quadratic  programming  problem  as 
will  be  shown  below. 

Method  of  Inversion 


Using  the  representation  theorem  for  seismic  sources,  a  connection  can  be  made  between 
source  parameters  and  observed  displacements  in  terms  of  equivalent  body  forces.  This  is 
a  force  system  which  would  produce  displacements  of  the  earth’s  surface  equivalent  to  that 
from  the  true  physical  situation  i.e.  a  heterogeneous  fault  system  with  a  non-linear  rheology. 
A  general  indigenous  source  can  be  expressed  in  terms  of  equivalent  body  forces  (Backus 
and  Mulcahv,  1976),  resulting  in  the  relationship 


uk(x,t)  =  I  f  Gki(x,t;r,t)fi(r,t)dVdt 

J  —  oo  JV 


where  uk  is  the  kth  component  of  displacement,  Gkl(x,t\r,t)  arc  the  Green  functions  con¬ 
taining  propagation  effects,  /,( r , t)  are  the  sum  of  the  equivalent  body  forces,  and  V  is 


the  source  region  where  the  /,(r,t)  are  non-zero.  The  summation  convention  is  assumed 
for  repeated  indices.  The  (7*,(x,  t:  r,  t)  may  be  expanded  in  a  Taylors  series  (Stump  and 
Johnson,  1977)  about  the  point  r  =  £, 


Gki(x,t-,r,t)  =  -^(n  -  ft) ..  .  (r„  -  £„)Gwj,...j„(r,  <;£,/) 

n! 

Define  the  force  moment  tensor 

=  JvSrji-ZJl)...(rjn-Z,n)ft(r,t)dV. 

and  the  displacement  becomes 
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n—  1 

where  X  represents  a  temporal  convolution. 

If  the  source  dimensions  are  small  compared  to  the  wavelengths  of  interest,  then  it  is 
necessary  to  keep  only  the  term  for  n  =  tin  the  previous  expansion  and  this  results  in 

ufc(x,<)  =  6’fci.J(x.f:0.0)  x  MtJ(0,t)  (1) 

for  £  =  0.  \\y  Fourier  transforming  this  convolution  a  matrix  equation  is  derived, 

«tfc(x,/)  =  (7Wj(x./;0,0)A/0(0./)  (2) 

for  each  frequency.  It  is  possible  to  solve  either  equation  ( 1 )  or  equation  (2)  for  the  moment 
tensor  in  the  time  or  frequency  domain  respectively.  In  order  to  solve  the  time  domain 

equation  (1)  one  may  assume  a  source  time  function  which  is  the  same  for  all  moment 

tensor  components,  converting  the  convolution  into  a  multiplication.  Alternatively,  the 
convolution  equation  may  be  written  as  a  matrix  equation  and  the  resulting  large  block 
Toplitz  system  of  equations  may  be  solved  by  a  least  square  error  algorithm  (Sipkin  1982). 

Another  approach,  the  one  taken  in  this  paper,  is  to  work  in  the  frequency  domain, 
solving  equation  (2).  Each  component  may  have  an  arbitrary  source-time  function  and  it  is 
only  necessary  to  solve  a  much  smaller  system  of  equations  successively  at  each  frequency  of 
interest.  This  system  of  equations  may  be  solved  by  least  squares,  minimizing  the  P  norm 
of  the  residuals, 

h 

nun[Y(uk  -  GkiMi)2} 


where  l  is  the  index  vector  (i,j).  The  total  number  of  station  components  is  denoted  by 
K .  The  index  1  indicates  north  while  indices  2  and  3  represent  east  and  down  respectively. 
The  vector  M  is  given  by  the  real  and  imaginary  components  of  the  moment  tensor: 

'  reA/n  N 
imMu 
reM2\ 
tmA/21 
re  A/22 
imA/22 
reM3l 
im  A/31 
reM32 
imM32 
re  A/33 
\  imM33 

This  results  in  the  generalized  inverse  solution  which  may  be  written  in  terms  of  the  singular¬ 
valued  decomposition.  In  addition  to  being  a  well-known,  efficient  procedure,  this  method 
allowrs  one  to  form  the  data  and  model  resolution  matrices  and  the  unit  covariance  matrix 
(Menke  1984). 

It  is  important  to  find  a  ’’best  fitting  model”  which  the  least  squares  solution  provides, 
but  it  is  also  useful  to  use  the  data  to  answer  more  specific  geophysical  questions.  For 
example,  is  a  pure  double-couple  source  sufficient  to  satisfy  the  data,  or  is  some  volume 
change  present  in  the  source?  To  answer  such  a  question  it  is  necessary  to  find  the  extremum 
of  a  linear  combination  of  the  moment  tensor  components,  the  sum  of  the  diagonal  elements 
of  the  moment.  Alternatively,  one  could  minimize  the  modulus  squared  of  this  quantity, 
the  approach  we  advocate.  The  search  for  the  extrema  is  subject  to  the  constraint  that  the 
data  Ufc  are  satisfied  within  some  confidence  intervals  c^.  Thus  we  seek 


min\AiMi) 

(3) 

or. 

min[MmT  A  mi  M{\ 

(4) 

subject  to 

Uk~  ffc  <  GktMi  <  uk  +  tk  k  =  l,K. 

(5) 

T  e  Ami  matrix  is  crucial  because  it  defines  the  properties  which  will  be  bounded.  By 
changing  the  elements  of  Ami  a  variety  of  different  source  characteristics  can  be  considered: 
the  largest  and  smallest  vertical  strike-slip  fault  source,  the  most  and  least  thrust-like  source, 
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etc.  Table  1  presents  the  coefficients  of  the  quadratic  form  ( Amt  in  equation  -1)  for  a  number 
of  source  properties.  By  changing  the  sign  of  the  elements  the  problem  is  changed  from  a 
minimization  to  a  maximization. 

The  system  of  equations  and  inequalities  given  above  results  from  the  local  imposition 
of  error  bounds  for  each  constraint  equation,  that  is,  the  confidence  interval  for  each  datum 
is  a  strict  interval  which  no  datum  may  exceed.  If  the  data  set  is  large  it  is  unlikely  that 
some  intervals  will  not  be  exceeded  (Oldenberg  1983).  As  an  alternative  to  this  imposition 
of  absolute  confidence  bounds  on  the  data  it  is  possible  to  impose  a  statistical  bound.  Using 
this  approach,  one  would  require  that  the  i  th  datum  is  satisfied  within  some  unspecified 
error  bound  e,.  Specifically,  the  sum  of  the  error  for  each  constraint  should  not  exceed  some 
pre  determined  value  E.  Algebraically,  instead  of  the  constraints  in  equation  (5)  we  have, 

Uk  <  OkiMi  +  ffc  k  =  1,  K 

CkiMi  -  ck  <  uk  k  =  1.  A'  (6) 

and 

*=l 

Ck  >  0  k  —  \.K 

,S  >  0 

The  additional  constraint  has  been  imposed  that  all  the  e,'s  and  s  are  always  positive.  If 
the  errors  e,  are  independent,  normal,  random  variables,  the  value  of  E  may  be  calculated 
using  a  priori  estimates  of  their  mean  and  standard  deviations  cr,  (Parker  and  McNutt 
1980).  Unfortunately,  the  statistical  parameters  of  errors  associated  with  waveforms  are 
not  known  a  priori  In  addition  to  microseismic  noise,  which  may  be  estimated  by  pre¬ 
event  noise  samples,  there  is  signal  generated  noise  due  to  errors  in  the  estimation  of  the 
Green  function.  This  noise,  in  both  phase  and  amplitude,  is  often  more  significant  than 
microseismic  noise  and  must  be  accounted  for.  We  will  estimate  these  bounds  by  taking  the 
difference  between  the  observed  data  and  data  predicted  by  a  model  derived  by  an  inversion 
procedure  such  as  least  squares.  Another  method  to  estimate  errors  in  the  waveforms  is 
through  Monte  Carlo  simulation  of  data  sets.  This  involves  using  the  least  squares  solution 
to  generate  a  data  set  which  is  then  perturbed  by  an  a  priori  distribution  of  errors.  The 
statistic  of  this  data  set  may  be  computed  and  used  to  compute  confidence  intervals  with 
which  to  derive  extreme  models.  Unfortunately,  present  day  models  of  lateral  heterogeneity 
are  not  adaquate  to  fully  represent  the  Green  function  errors. 

The  algorithm  used  to  minimize  equation  (3)  subject  to  the  constraints  of  equation 
(5)  or  (6)  is  the  well  known  simplex  algorithm  for  linear  programming.  To  minimize  the 
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Table  1.  Coefficients  for  the  quadratic  form  in  equation  (4)  for  a  variety  of  moment  tensor 
properties.  The  real  and  imaginary  components  have  identical  coefficients 
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quadratic  functional  in  equation  (4)  subject  to  the  constraints  requires  quadratic  program¬ 
ming  (Dantzig  1963).  Bv  finding  the  minimum  and  maximum  of  the  functional  subject 
to  the  constraints  it  is  possible  to  derive  bounds  on  model  properties  from  bounds  on  the 
data.  This  is  one  method  to  explore  the  range  of  possible  models  rather  than  derive  a  single 
solution. 

As  mentioned  previously,  formulating  the  problem  in  terms  of  the  square  of  the  modulus 
of  a  linear  sum  of  moment  tensor  components  allows  the  transition  between  moment  tensor 
properties  in  the  time  and  frequency  domains.  This  follows  from  Parseval’s  theorem  which 
states  that  the  integral  of  the  square  of  the  modulus  of  a  function  is  equal  to  the  integral  of 
the  square  of  the  modulus  of  the  Fourier  transform  of  the  function  (Bracewell  1978).  If  one 
considers  the  functional  F(«v, )  =  as  a  finite  Fourier  transform  of  the  time  series  /(f), 

then  the  square  of  the  modulus  of  F(uy)  may  be  represented  as  a  quadratic  form  involving 
a  vector  M  composed  of  the  real  and  imaginary  parts  of  the  moment  tensor  components. 

F(u.-,)FUr  -  MT AM.  (7) 

Where  F(w’,)*  is  the  complex  conjugate  of  F(~,).  Because  of  Parseval's  theorem,  the  sum 
of  the  function  in  equation  (7)  over  all  frequencies  has  the  same  value  as  the  sum  of  the 
amplitude  squared  of  the  time  series  /(f)  over  all  time.  Fortunately,  our  time  series  are  finite 
and  our  seismic  spectra  are  band-limited,  due  to  instrument  transfer  functions.  Bence,  the 
extremum  of  this  sum  is  the  same  in  the  time  and  frequency  domain.  Therefore,  it  is 
possible  to  compute  the  minimum  or  maximum  value  of  the  square  of  the  modulus  of  any 
linear  combination  of  moment  tensor  components  in  the  frequency  domain,  sum  the  extreme 
values  over  all  frequencies,  and  interpret  this  directly  in  the  time  domain.  In  what  follows 
we  will  consider  the  moment  rate  tensor  rather  than  the  moment  tensor.  The  rate  tensor 
has  no  static  offset  and  hence  returns  to  zero  over  time. 

Deep  Earthquakes 

IX  bate  has  continued  about  the  presence  of  an  isotropic  moment  tensor  component  in 
deep  and  intermediate  subduction  zone  earthquakes  (I)ziewonski  and  Gilbert  1974,  Okal 
and  Geller  1979.  Silver  and  Jordan  1982,  Hodder  1984.  and  Riedesel  and  Jordan  1985).  It 
seems  reasonable  to  expect  subducting  regions  to  be  areas  of  compaction  accompanied  by 
dewatering,  high  fluid  pressure,  crack  closure  and  phase  changes.  The  essential  question  is  if 
any  of  these  phase  changes  are  meta-stable.  Only  then  can  such  reactions  produce  rapid  fail¬ 
ure  which  excites  high-frequency  waves  in  the  Earth.  Recently,  high  pressure  experiments 
(Kirby  1987,  Meade  and  Jeanloz  1988)  have  detected  shear  instabilities  associated  with 
phase  changes.  The  most  recent,  work  recorded  sudden  failure  and  acoustic  emissions  when 
simulated  ocean  lithosphere  was  subjected  to  pressures  of  up  to  20  GPa  at  temperatures  be¬ 
low  900  deg  K .  If  such  physical  processes  are  occuring  t  hey  should  be  seismically  detectable. 
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Some  investigators  have  emphasized  the  possible  trade-off  between  lateral  heterogeneities 
and  a  possible  isotropic  component  (Okal  and  Geller  1979).  In  addition,  source  complexity 
such  as  curved  faults  and  multiple  rupture  can  give  rise  to  non  double-couple  moment  tensor 
solutions  (Sipkin  1986).  The  derivation  of  extremal  models  for  the  time-varying  moment 
tensor  seems  well  suited  to  address  the  question  of  the  presence  of  an  isotropic  moment 
component,  given  microseismic  noise,  complex  sources,  and  lateral  heterogeneities. 

With  this  in  mind  we  consider  a  467.7  km  deep,  magnitude  5.6  earthquake  near  the 
Bonin  Islands  recorded  by  15  GDSN  long  period  instruments  (Figure  1  shows  the  station 
distribution).  The  fifteen  stations  span  a  distance  range  from  38.5°  to  91.7°  but  have  a  gap 
in  coverage  in  the  southeast  quadrant.  Shown  in  Figure  2  are  the  15  long  period,  vertical 
GDSN  seismograms.  The  first  256  seconds  of  the  records  were  used  in  the  inversions.  The 
signal  to  noise  ratio  is  quite  high  (consider  the  pre-event  microseismic  noise)  and  the  arrivals 
coherent  across  the  array  of  stations. 

Tests  with  synthetic  data 

To  illustrate  the  method  and  also  explore  its  capabilities,  we  first  consider  some  tests 
with  synthetic  data.  We  adopt  the  same  number,  distribution,  and  type  of  stations  as 
for  the  Bonin  Islands  earthquake  (Figure  1),  but  assume  that  the  event  is  at  a  depth  of 
467.7  km  with  the  source  specified  by  the  time- varying  moment  rate  tensor  shown  in  Figure 
3.  This  is  a  dip-slip  event  with  a  superimposed  isotropic  component.  The  body  waves 
from  the  event  (P,  pP,  and  sP)  were  calculated  by  the  WKBJ  method  (Chapman,  1978) 
for  a  PREM  Earth  model  (Dziewonski  and  Anderson,  1981)  and  then  low-pass  filtered  to 
obtain  the  vertical-component  synthetic  seismograms  shown  in  Figure  4.  Random  noise 
with  an  amplitude  10%  that  of  the  maximum  signal  has  been  added  to  each  seismogram 
to  simulate  microseismic  noise  with  a  signai-to-noise  ratio  of  10:1.  This  example  is  only 
intended  to  simulate  microseismic  noise  which  should  be  uncorrelated  for  the  global  station 
distribution  considered.  This  is  not  the  case  with  signal  generated  noise  arising  from  lateral 
heterogeneities  which  may  be  correlated  when  receiver  crustal  structures  are  similar. 

We  now  ask  the  question,  what  can  be  said  about  the  possibility  of  an  isotropic  com¬ 
ponent  in  the  source  given  the  noisy,  band-limited  seismograms?  One  way  to  answer  this 
question  is  to  compute  the  least-squares  solution  at  each  frequency  and  then  Fourier  trans¬ 
form  the  results  into  the  time  domain  where  the  moment  tensor  trace  can  be  computed. 
Such  results  are  shown  in  Figure  5  where  it  can  be  seen  that  they  do  a  good  job  of  recover¬ 
ing  the  known  solution  of  Figure  3.  But  the  question  remains,  because  of  the  microseismic 
noise  in  the  seismogram,  what  is  the  range  of  possible  moment  tensor  traces?  One  way 
to  quantify  this  is  to  find  the  models  with  the  maximum  and  minimum  trace  amplitudes 
squared  and  still  satisfying  the  data  within  the  errors  introduced  by  the  noise.  That  is. 
minimize  the  functional  in  equation  (4)  with  the  elements  of  the  Ami  matrix  given  by  the 
explosive  terms  in  Table  1. 


9 


BONIN  ISLANDS  EVENT 
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Figure  1.  Station  distribution  of  GDSN  instruments  recording  the  Bonin  Islands  earthquake 
of  October  4,  1985.  The  earthquake  epicenter  is  denoted  by  a  cross  at  the  center  of  the  map 
and  the  stations  by  triangles.  The  magnitude  5.6  event  was  at  a  depth  of  467.7  km. 


Figure  2.  GDSN  long  period  seismograms  from  the  Bonin  Islands  event.  The  early  P  arrivals 
(P,  pP,  and  PP)  are  present  in  this  section. 
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Figure  3.  Moment  tensor  source  used  for  the  synthetic  modelling.  The  source  consists  of  a  dip- 
slip  event  with  a  superimposed  isotropic  component.  The  isotropic  component  constituted 
60  %  of  the  source. 
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Figure  4.  Synthetic  seismograms,  with  10  %  noise,  from  the  source  model  in  Figure  3  and 
recorded  at  the  stations  in  Figure  1.  The  VVKBJ  method  was  used  to  compute  the  seismo¬ 
grams  using  the  PREM  velocity  structure. 
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Figure  5.  Least  squares  solution  for  the  moment  rate  tensor  using  the  noisy  seismograms. 
The  source  moment  rate  tensor  from  Figure  3  has  essentially  been  recovered.  The  diagonal 
elements  and  A/33  are  superimposed  on  one  another. 
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To  use  the  quadratic  programming  methods  described  above  it  is  necessary  to  first 
estimate  bounds  on  errors  in  the  data.  In  this  example,  absolute  confidence  bounds  will 
be  used,  i.e.  constraints  (5)  will  be  considered.  The  bounds  were  computed  by  taking  the 
least-squares  solution  in  Figure  5  and  predicting  the  displacement  observed  at  each  station. 
The  difference  between  the  vector  of  observed  components,  u,  and  the  vector  of  predicted 
components,  up,  is  used  as  an  estimate  of  the  errors  for  the  data, 

(  =  U  -  Up 

In  analogy  with  the  95%  confidence  interval  of  normal  distributions,  twice  c  was  used  as 
confidence  intervals  on  the  data.  These  bounds  on  the  real  and  imaginary  components  of 
the  data  as  a  function  of  frequency  are  shown  in  Figure  6  for  the  third  station.  Using 
these  estimates  of  the  error  bounds  it  is  possible  to  compute  the  model  with  the  minimum 
squared  moment  tensor  trace  modulus  and  the  results  are  shown  in  Figure  7.  Since  this 
is  an  extreme  solution,  it  exhibits  considerably  more  deviation  from  the  known  solution 
than  does  the  least-squares  solution  of  Figure  5.  A  similar  computation  can  be  performed 
for  the  maximum  trace  modulus,  and  the  moment  rate  tensor  traces  for  the  two  extreme 
models,  the  maximum  and  and  minimum  sum  of  squared  traces,  are  compared  with  the 
least-squares  inverse  moment  rate  tensor  trace  in  Figure  8.  As  expected,  the  least  squares 
solution  lies  between  the  two  extreme  solutions.  What  the  extreme  solutions  provide  is  a 
measure  of  the  uncertainty  in  the  least-squares  solution  due  to  random  noise  in  the  data. 
The  effect  of  increasing  the  noise  in  the  data  is  shown  in  Figure  9.  Here,  the  minimum 
extreme  traces  are  shown  for  varying  signal  to  noise  ratios.  The  bounds  get  progressively 
wider  as  the  additive  noise  increases,  the  lower  bound  tending  toward  zero  with  greater 
noise.  This  shows  how  widening  confidence  bounds  on  the  data,  caused  by  a  decreasing 
signal  to  noise  ratio,  leads  to  wider  bounds  on  the  isotropic  moment  tensor  component. 

Lateral  heterogeneity  in  both  attenuation  and  velocity  structure  also  effects  moment 
tensor  estimation  because  it  introduces  phase  and  amplitude  errors  into  the  data.  If  the 
Green  function  used  differs  from  the  real  Earth,  differences  between  the  observed  data 
and  the  predicted  data  will  occur.  The  errors  in  the  Green  function  enter  the  above  data 
confidence  intervals  through  the  assumed  data  resolution  matrix  Ra-  This  is  a  matrix  which 
relates  the  observed  data  to  the  data  predicted  assuming  a  particular  velocity  structure 
(Green  furction)  such  as  PREM.  For  an  over-determined  problem  it  can  be  written  directly 
in  terms  of  the  assumed  Green  function  Ga  (Aki  and  Richards  1980,  Menke  1984), 

Ra  =  Ga  (Gj  Ga  )~'gJ  .  (8) 

Contained  in  Ra  are  the  effects  of  the  experimental  setup  (station  distribution,  microseismic 
noise  etc.)  as  well  as  the  effects  of  lateral  heterogeneity.  To  make  this  clearer  Ra  can  be 
written  in  terms  of  the  "true’'  Green  function  Gt  .  We  make  the  assumption  that  the  "true" 
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Figure  8.  Error  bounds  on  the  real  and  imaginary  components  of  the  spectra  for  the  record 
from  station  3.  Each  point  in  the  figure  is  at  a  different  frequency. 
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Figure  7.  Model  having  the  minimum  squared  modulus  for  the  moment  rate  tensor  trace. 
Essential  elements  of  the  original  source  model  axe  recovered  but  the  diagonal  elements  are 
no  longer  identical.  Furthermore,  the  A/21  and  the  M 3t  elements  are  no  longer  zero. 


17 


MOMENT  RATE  (XIOl 19  DYNE- CM/SEC) 


BONIN  TRACES  (0.1  NOISE) 


Figure  8.  Comparison  of  the  moment  rate  tensor  traces  of  the  three  models;  least  squares, 
maximum  total  squared  modulus,  minimum  total  squared  modulus. 
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Figure  9.  Ettect  oa  the  minimum  total  trace  squared  solution  of  varying  the  signal  to  noise 
ratio  of  the  seismograms.  It  is  seen  here  that  as  the  percentage  noise  increases  the  minimum 
approaches  zero. 
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Green  function  is  a  perturbation  of  the  assumed  Green  function, 


Gt  =  Ga  +7-  (9) 

We  assume  that  ||7||  <  ||Ga  ||.  the  matrix  norm  of  the  Green  function  perturbation  is 
much  less  than  the  matrix  norm  of  assumed  Green  function.  Substituting  equation  (9)  into 
equation  (8)  results  in  the  following  expression  for  Ra  , 

Ra  =  (Gt  -7)(G?  Gt  -GT7-7TGt  -7T7)_1(G?  -7T).  (10) 

Neglecting  terms  higher  than  first  order  in  7,  factoring  out  G^  Gt  and  expanding  the 
inverse  reoults  in, 

Ra  =  Gt  (G?  Gt  r'GT  +  Gt  (G?  Gt  )-2G?  7G?  +  Gt  (G?  Gt  )-27TGt  G? 
-7(G?  Gt  r1  Gj  -  Gt  (G?  Gt  )-‘7T 

The  first  term  on  the  right  is  the  resolution  matrix  in  terms  of  the  "true”  Earth  model  and 
Green  function,  Gt  .  Therefore,  the  confidence  bounds,  e  in  equation  (5)  are  given  by, 

<  =  (I  ~  Rt  -  Gt  (Gj  Gt  )"2GtT  7GJ  -  Gt  (GtT  Gt  r27TGt  Gj  (11) 

+7(G?  Gt  )-lG?  +Gt  (G?  Gt  )-,7T)u. 

Rt  does  not  depend  on  the  perturbations  of  the  Green  function,  only  on  the  structure  of 
the  experiment.  It  is  now  clear  that,  even  if  the  experiment  were  structured  such  that  the 
data  were  perfectly  resolved,  the  presence  of  lateral  heterogeneity,  nonzero  7,  could  still 
produce  non-zero  confidence  bounds  on  the  data.  Therefore,  if  the  differences  between  the 
observed  and  predicted  data  are  used  as  confidence  bounds  on  the  data,  then  the  effect  of 
model  errors  will  be  incorporated  in  the  extreme  model  estimates.  In  the  presence  of  lateral 
heterogeneity  the  extreme  bounds  on  the  model  properties  will  be  wider. 

A  synthetic  test  was  also  performed  to  explore  the  effect  of  the  lateral  heterogeneity.  The 
test  was  similar  to  the  previous  one  with  the  earthquake  at  625  km  depth  and  the  station 
distribution  as  shown  in  Figure  1.  The  mechanism  is  the  same  as  before  (Figure  3)  but  now 
the  earth  model  is  not  homogeneous.  Instead,  each  path  consists  of  a  perturbed  version  of 
PREM.  Specifically,  a  2%  random  perturbation  in  velocity  was  added  at  each  depth  in  the 
model  with  independent  perturbations  for  each  ray  path.  The  lateral  heterogeneity  results 
in  the  synthetic  seismograms  shown  in  Figure  10  with  errors  in  the  phase  and  amplitude 
present.  When  the  extremal  solutions,  minimum  and  maximum  trace  squared,  are  computed 
and  presented  with  the  least-squares  solution,  they  differ  substantially  (Figure  11).  By 
comparing  these  models  with  similar  models  derived  with  10%  random  noise  present  (Figure 
8)  it  is  seen  that  the  effect  of  lateral  heterogeneity  can  be  quite  strong. 


Figure  10.  Seismograms  with  perturbed  arrivals  due  to  lateral  velocity  inhomogeneities.  Two 
percent  random  perturbations  on  the  PREM  model  were  applied  to  each  source-receiver 
path 
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*  igure  11.  Moment  rate  tensor  traces  of  the  least  squares  and  maximum  and  minimum  ex¬ 
treme  trace  solutions.  These  solutions  are  the  result  of  inverting  the  perturbed  seismograms 
in  Figure  10  while  assuming  a  laterally  homogeneous  velocity  structure  (PREM). 
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These  examples  demonstrate  some  strengths  and  weaknesses  of  the  approach  taken 
so  far.  First,  the  error  bounds  derived  depend  directly  on  the  knowledge  of  the  velocity 
structure,  the  point  source  assumption,  and  the  microseismic  noise.  This  dependence  is 
through  the  least  squares  solution  which  is  used  to  construct  the  predicted  data.  For  an 
over-determined  problem,  as  our  knowledge  improves,  for  example  through  the  modeling 
of  lateral  heterogeneities,  the  bounds  on  the  data  will  narrow.  However,  with  the  absolute 
bounds  used,  the  strict  inequalities  of  equation  (5),  can  be  either  too  restrictive  or  too  wide. 
This  is  because  the  confidence  bounds  must  be  satisfied  exactly,  no  single  point  may  exceed 
the  bounds.  Therefore,  the  bounds  must  be  made  wide  in  order  for  the  probability  that  any 
data  exceeds  these  bounds  to  be  small  (Oldenburg  1983).  This  can  be  improved  by  the  use 
the  statistical  bounds  given  in  equation  (6). 

Results  for  the  Bonin  Islands  earthquake 

Given  the  results  of  these  synthetic  tests  which  illustrate  the  effects  of  random  noise  and 
Earth  heterogeneity,  we  can  now  return  to  the  inversion  of  the  data  shown  in  Figure  2  for 
the  Bonin  Islands  earthquake.  Recall  that  the  basic  question  to  be  answered  is  whether  the 
source  of  this  earthquake  had  a  significant  isotropic  component.  As  a  prelude  to  the  inversion 
of  the  actual  data,  one  further  set  of  synthetic  tests  was  performed.  Using  the  station 
distribution  and  event  location  identical  to  the  Bonin  Islands  earthquake  and  assuming  10% 
random  noise  and  2%  lateral  velocity  perturbations,  the  mimimum  trace  squared  solution 
was  obtained  for  a  series  of  sources  having  different  relative  sizes  of  the  isotropic  component. 
We  consider  an  isotropic  component  to  be  identifiable  if  it  is  distinguishable  from  other 
features  of  the  solution  which  are  caused  by  the  mapping  of  microseismic  noise  and  Green 
function  errors  into  the  solution.  The  results,  which  are  shown  in  Figure  12,  allow  one  to 
ask  the  question:  What  proportion  of  isotropic  component  must  be  present  in  the  source 
in  order  for  it  to  be  unambiguously  identified  in  the  inversion  results?  The  answer  is  fairly 
high:  more  than  10%>  of  the  source  has  to  be  due  to  volume  change  alone.  When  the  same 
exercise  was  conducted  using  perturbations  of  5%  at  least  20%  of  the  source  had  to  be 
isotropic  for  detection.  Thus  it  may  not  be  possible  to  discern  a  small  isotropic  moment 
tensor  without  more  and  better  data  and  without  better  modelling  of  the  velocity  structure. 

Turning  now  to  the  data  shown  in  Figure  2,  the  least  squares  solution  is  shown  in 
Figure  13.  The  main  part  of  the  source  is  an  initial  pulse  of  about  20  sec  duration,  which 
reflects  the  bandwidth  of  the  instrumentation.  For  this  initial  pulse,  the  principal  axes 
of  the  deviatoric  part  of  the  source  have  the  approximate  orientations  (plunge,  azimuth): 
tension  axis  =  (10,  42);  intermediate  axis  =  (45,  144);  compression  axis  =  (40,  320).  This 
is  similar  to  the  Harvard  solution  (Bull.  ISC):  tension  axis  =  (16,  49);  intermediate  axis 
=  (23,  146);  compression  axis  =  (61,  287).  It  is  also  generally  consistent  with  the  fault 
plane  solutions  of  other  earthquakes  in  this  region  (Burbach  and  Frolich,  1986).  However, 
the  primary  interest  in  this  study  is  the  trace  of  the  moment  tensor  and  it  is  obvious  that 
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Figure  12.  Simulation  of  the  moment  rate  traces  in  the  minimum  total  trace  squared  solu¬ 
tion  for  varying  proportions  of  the  isotropic  component.  The  velocity  structure  was  again 
assumed  to  be  PRJEM  while  the  structure  used  to  generate  the  seismograms  contained  two 
percent  perturbations  of  PREM.  The  proportion  of  isotropic  component  is  given  by  the  ratio 
of  the  isotropic  component  to  the  sum  of  the  absolute  values  of  the  moment  components. 
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this  is  significantly  different  from  zero  for  the  least-squares  solution,  because  in  Figure  13 
the  M 33  component  is  much  larger  than  the  Mu  and  M 22  components.  This  trace  from 
the  least-sqaures  solution  is  shown  in  Figure  14  along  with  maximum  and  minimum  trace 
solutions  obtained  with  the  quadratic  programming  approach.  The  error  bound  assumed  in 
the  extremal  solutions  was  twice  the  root -mean- square  error  of  the  least-squares  solution. 
The  least-squares  solution  has  a  prominent  negative  isotropic  component  in  the  first  20  sec 
which  could  be  interpreted  as  a  volume  change  at  the  source.  Unfortunatelj’,  the  mimimum 
extreme  solution  shows  that,  given  the  uncertainty  in  the  data  and  the  Green  functions 
used  in  the  inversion,  this  volume  change  suggested  by  the  least-squares  solution  may  not 
be  significant.  On  the  basis  of  the  synthetic  tests  and  the  fairly  large  signal  to  noise  ratio 
in  the  data  (Figure  2).  it  seems  likely  that  lateral  heterogeneity  is  largely  responsible  for 
the  width  of  the  bounds  and  therefore  the  weakness  of  any  conclusion  which  may  be  drawn 
from  these  results. 

In  order  to  test  if  the  minimum  volume  change  is  significantly  different  from  zero,  the 
best  fitting  solution  with  the  trace  constrained  to  be  zero  could  be  found.  Then  a  Monte 
Carlo  simulation,  considering  all  error  sources,  could  be  used  to  compile  a  population  of 
waveform  data  sets.  An  inversion  of  these  data  sets  would  give  statistics  on  the  zero  trace 
solutions.  This  would  allow  one  to  statistically  test  if  the  minimum  squared  trace  solution  is 
significantly  different  from  zero.  We  feel  that,  at  present,  models  of  velocity  and  attenuation 
lateral  heterogeneity  are  not  yet  adaquate  for  this.  Instead,  we  rely  on  a  comparison  of  the 
isotropic  component  with  the  other  moment  tensor  components  to  estimate  significance. 

Nuclear  Explosion  Sources 

One  seismic  source  which  is  known  to  have  a  large  isotropic  component  is  a  nuclear 
explosion.  It  is  an  ideal  case  for  the  computation  of  extremal  moment  tensor  models.  Near 
source  data  from  a  nuclear  explosion  can  be  used  to  compute  upper  and  lower  hounds  on 
squared  moment  tensor  trace,  which  is  a  measure  of  the  volume  change  associated  with  the 
source.  Thus  the  solutions  are.  respectively,  the  most  and  least  explosion-like  solutions.  This 
has  important  applications  in  the  verification  of  nuclear  explosions  because  these  solutions 
provide  best  and  worst  cases  by  which  to  decide  if  an  event  was  a  nuclear  explosion. 

The  Harzer  experiment  of  June  6,  1981,  a  nuclear  explosion  of  equivalent  magnitude  5.5 
at  a  depth  of  637  m  was  recorded  by  eight,  three  component,  broadband,  digital,  accelerom¬ 
eters.  The  details  of  the  collection  and  interpretation  of  the  data  can  be  found  in  Johnson 
(1988)  and  will  only  be  briefly  reviewed  here.  The  near  source  network  was  azimuthally 
distributed  around  the  epicenter  with  stations  from  2.4  km  to  6.6  km  from  the  event  (Figure 
15).  The  velocity  records  are  shown  in  Figure  16.  The  network  was  located  in  the  Silent 
Canyon  Caldera,  a  heterogeneous  velocity  structure.  The  averaged  one-dimensional  velocity 
structure  in  this  Caldera  has  been  studied  by  Leonard  and  Johnson  (1987)  and  their  model 
was  used  in  moment  tensor  inversions.  A  modified  reflectivity  method  (Kind  1978)  was 
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Figure  15.  Harzer  experimental  setup.  Station  distribution  within  the  Silent  Canyon  Caldera. 
The  stations  were  three  component  broadband  accelerometers. 
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Figure  10.  Rotated  velocity  records  from  the  Harzer  explosion.  Note  the  significant  energy 
on  the  transverse  components. 
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used  to  compute  the  Green  functions. 

Because  of  difficulties  at  one  of  the  sites,  seven  stations  were  used  in  the  inversion,  giv¬ 
ing  a  total  of  21  components.  The  least-squares  solution  for  the  moment  rate  tensor  of  this 
overdetermined  problem  is  shown  in  Figure  17.  The  most  important  part  of  the  moment  rate 
tensor  is  the  initial  short-period  pulse  which  is  followed  by  longer-period  oscillations  which 
are  less  coherent  and  more  poorly  resolved.  This  initial  pulse  is  most  pronounced  on  the 
trace  components  Mu,  A/221  and  A/33,  although  there  is  still  energy  on  the  off-diagonal  ele¬ 
ments,  such  as  A/23.  The  A/n  and  A/22  elements  are  fairly  similar  in  their  time  dependence, 
but  the  A/33  element  is  somewhat  different,  containing  a  large  secondary  pulse  at  about  2 
seconds.  When  this  solution  is  combined  with  the  Green  functions,  predicted  seismograms 
are  obtained  which  do  a  reasonably  good  job  of  explaining  the  major  features  of  the  ob¬ 
served  data  in  Figure  16,  with  average  correlation  coefficients  of  about  0.5.  However,  there 
still  remain  significant  differences  between  the  predicted  and  observed  seismograms,  partic¬ 
ularly  on  the  transverse  components.  Much  of  this  difference  can  probably  be  attributed  to 
effects  which  were  not  taken  into  account  in  the  modelling,  such  as  lateral  heterogeneities 
and  scattering  in  the  wave  propagation  and  secondary  source  effects  such  as  spall. 

Given  that  the  source  has  been  less  than  perfectly  resolved  by  the  least-squares  solution 
for  the  moment  tensor,  what  can  be  said  about  the  uncertainty  in  the  explosive  part  of 
the  source?  This  question  is  answered  in  Figure  18  which  shows  the  least-squares  solution 
for  the  moment  rate  ten.  or  trace  along  with  the  estimates  for  the  minimum  and  maximum 
moment  rate  tensor  traces.  It  is  clear  that  an  initial  compressional  pulse  is  present  on 
all  three  solutions.  Thus,  even  in  the  presence  of  fairly  large  error  bounds  due  to  both 
random  noise  and  deficiencies  in  the  modelling  process,  it  can  be  concluded  that  this  source 
has  a  clear  explosive  component.  The  extremal  solutions  are  also  useful  in  associating  an 
uncertainty  with  the  first  pulse  on  the  trace,  which  is  directly  related  to  the  yield  of  the 
explosion. 

Up  to  this  point  all  of  the  calculations  have  employed  the  local  error  bounds  of  equa¬ 
tion  (5),  but  there  are  situations  where  the  global  error  bounds  of  equation  (6)  might  be 
preferred.  The  results  of  using  these  two  types  of  error  bounds  are  compared  in  Figure  19 
for  the  minimum  trace  solution.  In  this  particular  case  the  choice  of  bound  does  not  cause 
enough  difference  in  the  results  to  affect  any  of  the  conclusions  based  upon  them. 

Conclusions 

A  method  has  been  developed  by  which  extreme  models  of  the  time-varying  moment 
tensor  may  be  constructed.  The  method  has  the  potential  to  answer  many  interesting 
geophysical  questions.  In  particular,  it  is  now  possible  to  assess  the  presence  or  absence  of 
volume  change  (isotropic  component)  in  seismic  sources.  It  has  been  comonly  assumed  that 
the  isotropic  component  of  the  moment  tensor  vanishes.  We  believe  that  this  assumption 
has  not  been  adaquately  examined  and  that  this  method  can  be  used  for  this  purpose. 
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Figure  19.  Minimum  moment  rate  trace  squared  solutions  computed  using  the  strict  data 
bounds  of  equation  (5)  and  the  global  data  bounds  of  equation  (6).  Again,  both  models 
have  an  obvious  isotropic  component. 


The  two  applications  described  above  illustrate  two  areas  in  which  extreme  bounds  on  the 
isotropic  component  are  useful:  deep  earthquakes  and  nuclear  explosions.  Surely,  many 
others  are  possible.  One  exciting  aspect  of  the  technique  is  that  it  makes  full  use  of  the 
waveforms.  With  the  advent  of  new  networks  of  wide  dynamic  range,  broad  band  digital 
seismometers  such  as  GEOSCOPE  and  IRIS  greater  resolution  will  be  possible.  As  can 
be  seen  from  the  Harzer  inversions,  even  high  frequency  data  can  give  exellent  inversion 
results. 

The  method  described  in  this  paper  provides  an  extreme  estimate  of  some  aspect  of 
the  moment  tensor,  given  the  recorded  seismograms  and  an  estimate  of  their  uncertainty  in 
either  the  time  domain  or  frequency  domain.  This  uncertainty  should  include  all  possible 
sources  of  error,  including  random  noise,  the  source  location,  the  earth  model  used  in 
calculating  the  Green  functions,  and  the  method  used  to  calculate  the  Green  functions.  In 
some  situations  it  may  be  possible  to  make  a  priori  estimates  of  all  these  errors,  but  in 
general  this  will  not  be  possible.  In  this  paper  we  have  proposed  the  alternative  procedure 
of  using  the  residual  of  the  least -squares  solution  as  an  estimate  of  the  total  error.  While 
this  procedure  has  the  undesireable  feature  that  the  error  estimate  depends  upon  the  data, 
it  does  provide  a  procedure  that  can  be  used  in  all  instances  and  it  seems  to  have  given 
reasonable  results  in  the  two  cases  considered.  The  method  is  flexible  in  the  sense  that  the 
error  constraints  can  be  applied  in  either  a  local  or  global  sense. 
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ABSTRACT 

The  Geysers  geothermal  field  is  the  site  of  intense  microseismicity  which 
appears  to  be  associated  with  steam  production.  It  appears  that  focal  mechanisms  of 
earthquakes  at  The  Geysers  vary  systematically  with  depth,  but  P-wave  first-motion 
focal  mechanism  studies  have  been  hampered  by  inadequate  resolution.  In  this  study 
an  unconstrained  frequency  domain  moment  tensor  inversion  method  is  used  to  over¬ 
come  P-wave  first-motion  focal  sphere  distribution  problems  and  to  investigate 
microearthquake  source  properties  .  A  goal  was  to  investigate  the  feasibility  of  using 
waveforms  to  invert  for  the  second-order  moment  tensor  of  microearthquakes  in  the 
complex  setting  of  The  Geysers.  Derived  frequency-domain  moment  tensors  for  two 
earthquakes  were  verified  by  mechanisms  estimated  from  P-wave  first  motions  and 
required  far  fewer  stations.  For  one  event,  19  P-wave  first  motions  were  insufficient 
to  distinguish  between  normal-slip  and  strike-slip  focal  mechanisms,  but  a  well  con¬ 
strained  strike-slip  solution  was  obtained  from  the  waveform  principal  moment  inver¬ 
sion  using  data  from  6  stations.  Improved  waveform  focal  mechanism  resolution  was 
a  direct  consequence  of  using  P  and  S-wave  data  together  in  a  progressive  velocity- 
hypocenter  inversion  to  minimize  Green  function  errors.  The  effects  of  hypocenter 
mislocation  and  velocity  model  Green  function  errors  on  moment  tensor  estimates 
were  investigated.  Synthetic  tests  indicate  that  these  errors  can  introduce  spurious  iso- 
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tropic  and  CLVD  components  as  large  as  26%  for  these  events  whereas  principal 
moment  orientations  errors  were  <8°.  In  spite  of  unfavorable  recording  geometries 
and  large  (0.6  km)  station  elevation  differences,  the  results  indicate  that  waveform 
moment  tensor  estimates  for  microearthquake  sources  can  be  robust  and  constrain 
source  mechanisms  using  data  from  a  relatively  small  number  of  stations. 


1.  Introduction 

Estimates  of  seismic  source  properties  are  among  the  most  important  pieces  of  information 
extracted  from  recordings  of  microearthquakes.  Source  studies  using  microcarthquakes  can  constrain 
the  mechanism  and  geometry  of  faulting  and  provide  estimates  of  principal  stress  orientations.  These 
estimates  form  a  basis  for  interpreting  deformation  associated  with  the  earthquakes  and  possible  rela¬ 
tions  between  seismicity  and  tectonic  stresses. 

The  Geysers  is  the  world’s  largest  generator  of  electricity  using  geothermal  energy  and  the  site  of 
intense  microseismicity  The  rate  of  seismicity  at  The  Geysers  is  45  times  the  regional  rate  (Ludwin  et 
al.,  1982).  Most  earthquakes  occur  within  or  just  below  the  shallow  steam  production  zone  at  depths  of 
2  to  4  km  below  the  mean  surface  elevation.  At  The  Geysers  the  relationship  between  seismicity,  tec¬ 
tonic,  and  locally  induced  stresses  is  unclear.  Although  it  appears  that  seismicity  at  The  Geysers  is 
induced  by  steam  production  activities,  a  specific  mechanism  has  not  been  determined  (Oppenheimer, 
1986;  Eberhart-  Phillips  and  Oppenheimer,  1984;  Bufe  et  al.,  1981).  Bufe  et  al.  (1981)  suggested  that 
the  wide  range  in  fault  plane  solutions  found  al  The  Geysers  were  a  function  of  ume. 

Oppenheimer  (1986)  estimated  the  stress  field  orientation  at  The  Geysers  from  210  fault  plane 
solutions.  Based  on  agreement  of  the  extensional  principal  stress  direction  estimated  from  the  seismicity 
as  a  whole  with  that  obtained  from  regional  geodetic  data  (Prescott  and  Yu,  1986),  he  concluded  that 
regional  tectonic  stresses  are  much  larger  than  the  stresses  induced  locally  through  geothermal  activities. 
Oppenheimer  (1986)  found  that  shallow  earthquake  focal  mechanisms  are  dominantly  strike-slip  and 
reverse  whereas  deeper  focal  mechanisms  predominantly  exhibit  normal  faulting.  He  concluded  that 
focal  mechanisms  at  The  Geysers  are  a  function  of  depth.  His  estimation  of  stress-field  mentation  and 
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variation  of  focal  mechanisms  with  depth  were  hampered  by  the  common  problem  of  nonunique  fault 
plane  solutions. 

For  many  shallow  events  there  is  an  ambiguity  between  pure  strike-slip  and  pure  dip-slip  mechan¬ 
isms.  These  ambiguities  are  not  due  to  a  paucity  of  P-wave  first  motion  data.  Numerous  P-wave  first 
motions  are  available  for  earthquake  at  The  Geysers  because  an  extensive  seismic  recording  network 
has  been  operated  since  1976  in  The  Geysers  area  by  the  U.S.  Geological  Survey  (USGS).  Often  how¬ 
ever,  P-wave  first  motions  are  absent  from  the  central  portions  of  the  upper  focal  hemisphere  for  shal¬ 
low  events  because  no  stations  are  close  enough  to  the  epicenter.  This  problem  persists  despite  an  aver¬ 
age  station  spacing  of  about  2  km  in  the  source  area  and  can  be  attributed  primarily  to  two  factors. 
The  earthquakes  are  very  shallow,  generally  2  to  4  km  below  mean  station  elevations,  and  velocity  gra¬ 
dients  (Figure  (1))  are  large.  These  two  factors  combine  to  severely  reduce  the  number  of  observations 
in  the  central  portion  of  the  upper  focal  hemisphere.  Observations  from  distant  stations  that  sample  the 
central  portion  of  the  lower  focal  hemisphere  are  absent  due  to  attenuation  of  microearthquake  signals. 
Oppenheimer  (1986)  noted  that  for  some  events,  fault  plane  solutions  were  completely  ambiguous; 
strike-slip,  reverse-slip,  and  normal-slip  solutions  could  fit  the  same  first  motion  data. 

An  alternative  approach  to  estimate  source  mechanisms  and  principal  moment  orientations  is  to 
use  a  waveform  inversion  method  to  estimate  seismic  moment  tensors  of  microearthquakes  at  The 
Geysers.  The  frequency  domain  method  of  Stump  and  Johnson  (1977)  is  used  here  to  estimate  second- 
order  moment  tensors  for  three  microearthquakes  at  The  Geysers  geothermal  field.  A  goal  is  to  deter¬ 
mine  the  utility  of  waveform  moment  tensor  inversion  to  supplement  P-wave  first-motion  data,  when 
they  are  of  insufficient  quantity  and  distribution  to  constrain  microearthquake  focal  mechanisms.  The 
primary  utility  of  the  moment  tensor  inversion  approach  in  this  context  is  to  obtain  well  constrained 
focal  mechanism  estimates  using  full  waveform  data  from  a  limited  number  of  stations. 

We  also  wanted  to  determine  the  reliability  of  microearthquake  waveform  moment  tensor  inver¬ 
sion  results.  Saikia  and  Herrmann  (1986)  could  not  independently  verify  the  accuracy  of  waveform 
moment  tensor  inversion  for  Arkansas  microearthquakes  because  first-motion  data  were  only  available 
from  a  small  number  of  stations.  The  large  number  of  P-wave  first  motions  available  for 
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microearthquakes  at  The  Geysers  provide  sufficient  focal  mechanism  constraints  to  check  waveform 
moment  tensor  estimates  for  two  of  the  earthquakes. 

In  contrast  to  teleseismic  and  regional  geometries,  microearthquake  mislocation  can  produce  sub¬ 
stantial  errors  in  distances,  azimuths,  and  take-off  angles  used  to  calculate  Green  functions.  An 
appraisal  methodology  is  developed  to  quantify  the  effects  of  these  components  of  Green  function  errors 
and  velocity  model  errors  on  moment  tensor  decompositions.  We  begin  by  outlining  the  source  charac¬ 
terization  and  moment  tensor  inversion  method  used.  Next,  the  data  set  and  details  of  the  inversions 
are  presented.  The  moment  tensor  error  appraisal  method  is  then  presented  and  applied  to  the  results  of 
moment  tensor  inversions. 

2.  Source  Characterization 

The  moment  tensor  formulation  is  used  to  represent  the  seismic  source  in  space  and  time.  Utiliz¬ 
ing  the  fact  that  a  seismic  source  can  be  represented  as  a  set  of  equivalent  body  forces,  the  source  can 
be  expanded  as  a  series  of  moments.  For  small  sources  or  large  wavelengths,  only  the  first  term  of  the 
series  is  retained  (point  source  approximation),  and  the  displacement  at  any  point  and  time  can  be  writ¬ 
ten  as 

f/t(x',rO  =  Gb.y(x\/';0,o)  8M1/(0,t')  (1) 

where  Uk  is  the  displacement  in  the  k  direction,  Gb  is  the  Green  function,  Mt/  is  the  moment  tensor,  . 
indicates  differentiation  with  respect  to  xf,  and  9  represents  temporal  convolution.  A  more  complete 
derivation  of  (1)  is  given  in  Stump  and  Johnson  (1977). 

In  the  frequency  domain,  equation  (1)  reduces  to 

£/*(£'/)  =  GUj(x'S;0.0)M„(QJ)  (2) 

If  the  propagation  paths  effects  (Gbj)  are  known,  one  can  determine  the  source  (Af„ )  from  a  set  of 
observational  data  ({/*)  by  solving  this  set  of  linear  equations. 

In  the  implementation  used  here,  Fourier  transforms  of  the  data  and  Green  functions  are  calcu¬ 
lated,  and  the  moment  rate  tensor  (M,,)  is  solved  for  in  the  frequency  domain.  An  inverse  Fourier 
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transform  is  used  to  obtain  MtJ  in  the  time  domain.  Then  is  detrended  to  eliminate  spurious  dc 
offsets.  The  resulting  estimates  of  are  integrated  to  yield  Mi; .  Detrending  of  Af,,  is  physically 
justified  because  the  moment  rate  tensor  elements  cannot  have  a  permanent  dc  offset.  If  perfect  data 
were  available,  detrending  would  not  be  required  but  all  seismic  data  are  intrinsically  bandlimited,  and 
the  instruments  used  here  (4.5  Hz  velocity  transducers)  have  limited  low  frequency  responses.  Detrend¬ 
ing  as  applied,  is  simply  a  high-pass  filter  operation  to  remove  spurious  low  frequency  (<  1  Hz)  noise. 

Since  the  complex  frequency  dependence  is  obtained  for  each  moment  tensor  element,  moment 
tensor  elements  are  not  required  to  have  a  common  time  function.  This  allows  inversion  for  complex 
sources  that  could  have  several  physical  source  components  with  different  time  histories.  An  alternative 
approach  is  to  solve  for  the  moment  tensor  element  time  functions  using  the  multichannel  vector 
decomposition  method  developed  by  Oldenburg  (1982),  as  presented  by  Sipkin  (1982).  In  cases  with 
source  multiplicity,  allowing  all  moment  tensor  elements  to  have  their  own  time  functions  eliminates 
errors  in  moment  tensor  estimates  that  are  inherent  in  time  domain  approaches  that  assume  a  common 
time  function  for  all  moment  tensor  elements  (see  Sipkin  (1986)  for  some  examples). 

The  moment  tensor  characterization  of  seismic  sources  provides  a  means  for  estimating  source 
properties  of  microearthquakes.  Stump  and  Johnson’s  (1977)  approach  is  completely  general;  no  restric¬ 
tive  assumptions  are  required  about  physical  source  types  or  the  time  dependence  of  moment  tensor  ele¬ 
ments.  All  physical  source  types  can  be  analyzed  including  isotropic  (volume)  sources,  compensated 
linear  vector  dipole  (CLVD)  sources,  and  double-couple  sources.  Any  of  these  source  types  can  be 
investigated  using  graphical  first-motion  methods  by  relaxing  the  common  double-couple  source 
assumption.  However,  the  linear  programming  moment  tensor  inversion  method  of  Julian  (1986)  pro¬ 
vides  the  most  powerful  means  to  utilize  first-motion  data. 

The  primary  advantage  of  the  using  moment  tensor  waveform  inversion  in  conjunction  with  first 
motion  information  is  that  fewer  recording  stations  are  needed  to  constrain  seismic  source  properties. 
Only  6  components  of  ground  motion  (two  three-component  stations)  are  required  in  theory,  although  in 
practice  about  15  components  of  ground  motion  are  recommended  to  ensure  reliable  results.  Another 
advantage  is  that  azimuth  and  takeoff  angle  coverage  need  not  be  as  comprehensive  as  when  only  P- 
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wave  first-motion  data  are  utilized.  Consequently,  it  is  possible  to  constrain  source  properties  of  earth¬ 
quakes  that  are  not  completely  surrounded  by  recording  stations,  something  that  is  not  generally  possi¬ 
ble  when  using  only  P-wave  first  motion  data. 

Frequency  domain  moment  tensor  inversion  has  not  been  applied  to  microearthquakes  before. 
Stump  and  Johnson  (1984)  have  used  the  method  to  characterize  nuclear  explosion  sources  using  near¬ 
field  data.  Time  domain  moment  tensor  inversions  with  restrictions  on  physical  source  type  (pure  devia- 
toric)  and  moment  tensor  time  dependence  have  been  applied  to  microearthquakc  data  by  Saikia  and 
Herrmann  (1986).  However,  their  approach  requires  assuming  that  all  moment  tensor  elements  have  the 
same  time  function  and  that  the  time  function  is  known  or  can  be  estimated.  The  result  of  this  type  of 
inversion  is  a  static  estimate  of  the  moment  tensor  elements.  Since  the  source  time  function  is  intrinsi¬ 
cally  unknown,  any  errors  in  the  assumed  time  function  will  produce  errors  in  the  static  moment  tensor 
estimate.  Further,  if  all  moment  tensor  elements  do  not  actually  have  the  same  lime  function,  another 
component  of  error  will  be  added  to  the  static  moment  tensor  estimate.  The  approach  used  here  allows 
each  moment  tensor  to  have  an  independent  time  function.  This  requires  more  data  than  the  lime 
domain  approach  of  Saikia  and  Herrmann  (19 86),  Langston  (1981),  and  Langston  and  Helmberger 
(1977),  but  yields  more  complete  information  about  source  properties. 

Estimated  moment  tensors  can  be  decomposed  into  isotropic  and  deviatoric  components.  The  rela¬ 
tions  are 

M-nt  (isotropic)  =  y  Mu  (3) 

Dt/  (deviatoric )  =  Af,y  -  Mn  b,, .  (4) 

If  prior  knowledge  is  available  about  the  source,  then  appropriate  constraints  can  be  placed  on  equation 
(2).  Constraints  on  equation  (2)  were  not  used  when  inverting  for  microearthquake  moment  tensors. 
Non-double-couple  earthquake  mechanisms  cannot  be  excluded  at  The  Geysers.  Consequently,  we 
wanted  to  investigate  if  unconstrained  moment  tensor  inversions  would  produce  moment  tensor  esti¬ 
mates  consistent  with  the  common  assumption  for  earthquakes  of  a  single  double-couple  source.  The 
frequency  domain  approach  was  used  to  avoid  errors  due  to  possible  source  multiplicity.  The  point 
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source  assumption  is  valid  for  the  microearthquakes  used  here,  since  source  dimensions  are  small  com¬ 
pared  to  the  wavelengths  represented  in  the  observed  data. 

The  eigenvalues  and  corresponding  eigenvectors  of  DtJ  describe  the  magnitude  and  orientation, 
respectively,  of  the  principal  moment  axes  (neglecting  gravity)  acting  at  the  source  .  These  principal 
moment  axes  represent  the  quantity  that  is  uniquely  determined  (within  a  range  of  uncertainties  due  to 
errors  in  Uk  and  Gbj)  by  moment  tensor  inversion.  Decomposition  of  DkJ  into  physical  source  com¬ 
ponents  is  fundamentally  nonunique  (Geller,  1976).  Julian  (1986)  used  linear-programming  methods  to 
investigate  the  range  of  possible  physical  source  mechanisms  that  are  consistent  with  a  particular  first- 
motion  or  amplitude  data  set.  The  common  approach  of  decomposing  D,;  into  double  couple  and 
CLVD  components  is  not  particularly  meaningful  due  to  its  intrinsic  nonuniqueness  unless  it  is  believed 
that  both  components  are  truly  contained  in  the  seismic  source.  A  simple  shear  dislocation  earthquake 
source  can  have  nonzero  CLVD  components  if  the  rupturing  fault  plane  has  nonzero  curvature  and  a 
nonzero  isotropic  component  if  a  fault  surface  is  irregular  or  imbedded  in  an  anisotropic  material 
(Backus  and  Mulcahy,  1976).  The  decomposition  of  D,,  into  CLVD  and  double  couple  components 
does  give  a  measure  of  the  departure  of  the  estimated  source  from  a  planar  faulting  single  double¬ 
couple  earthquake  model.  A  simple  measure  of  the  departure  of  D,j  from  a  single  double  couple  is  the 
ratio  of  the  smallest  and  largest  eigenvalues  of  0,, . 

A  simply  relationship  does  not  exist  between  principal  moment  and  principal  stress  orientations. 
The  maximum  principal  stress  orientation  is  poorly  resolved  by  the  principal  moment  orientations 
(McKenzie,  1 969).  Principal  moment  estimates  from  many  earthquakes  can  be  incorporated  into  stress 
tensor  inversions  (Angelier,  1984;  Gephart  and  Forsyth,  1984;  and  Michael,  1987).  Oppenheimer 
(1986)  used  the  method  of  Angelier  (1984)  to  estimate  principal  stress  orientations  at  The  Geysers  using 
focal  mechanism  data.  Michael  (1987)  demonstrated  that  principal  moment  orientation  data  can  provide 
sufficient  information  to  find  the  best  stress  tensor,  albeit  at  decreased  resolution  compared  to  using 
slickcnsidc  or  known  fault  orientation  data. 


3.  Data  Analysis 


A  9  station  three-component  digital  seismometer  network  was  deployed  in  The  Geysers  from  July 
21  to  August  15,  1982.  The  network  was  located  entirely  inside  the  geothermal  production  zone  to 
avoid  the  effects  of  significant  lateral  velocity  variations  (Ebcrhart-Phillips,  1986)  outside  this  zone.  A 
complete  description  of  the  recording  network  and  data  can  be  found  in  O’Connell  (1986).  Data  were 
recorded  at  200.32  samples/sec  using  three-component  4.5  Hz  velocity-transducer  geophones.  Before 
sampling,  the  data  were  anti-alias  lowpass  filtered  using  a  5-  pole  Butterworth  filter  with  a  50  Hz  comer 
frequency  and  high  pass  filtered  with  two  1-pole  Butterworth  filters  with  0.2  Hz  comers.  The  rapid 
decrease  of  displacement  magnification  of  the  velocity  transducers  below  4.5  Hz,  combined  with  the 
12-bit  resolution  of  the  recording  system,  limited  the  usable  frequency  band  to  the  range  from  1.0  to  50 
Hz.  Also,  electrical  problems  with  some  of  the  recorders  resulted  in  increased  noise  at  low  (<1  Hz)  fre¬ 
quencies,  so  moment  tensor  estimates  below  1  Hz  are  unreliable.  This  does  not  significantly  effect  the 
results  of  the  moment  tensor  inversions  because  source  comer  frequencies  for  the  microearthquakes 
used  here  are  between  6  and  10  Hz.  It  does  however,  necessitate  detrending  of  the  moment  rate  tensor 
in  the  time  domain,  as  discussed  earlier. 

The  maximum  observed  comer  frequencies  at  The  Geysers  is  strongly  correlated  with  particular 
station  sites.  Certain  stations  consistently  produced  lower  comer  frequencies  and  larger  high  frequencies 
rolloffs,  independent  of  epicentral  distance  or  hypocentra!  depth.  This  strongly  suggested  that  an  /m„ 
effect  (Hanks,  1982)  was  controlling  the  maximum  observed  comer  frequencies  at  some  stations.  An 
/mi,  of  15  to  20  Hz  was  observed  for  station  sited  on  slope  failure  materials.  To  reduce  problems  asso¬ 
ciated  with  variations  of  /m„  as  a  function  of  station  site,  the  data  were  lowpass  filtered  using  a  2-pole 
Butterworth  filter  with  a  10  Hz  comer  frequency.  This  also  decreased  the  burden  of  Green  function 
computations  by  reducing  the  maximum  frequency  required.  The  price  paid  is  that  moment  tensor  time 
functions  represent  lowpass-filtered  versions  of  the  true  source  time  functions. 

Green  functions  were  calculated  using  a  spectral  wavenumber-frequency  approach  similar  to  the 
reflectivity  method  of  Kind  (1978).  The  entire  model  between  the  free  surface  and  the  model  bottom 
constitutes  the  reflectivity  zone;  all  reverberations,  including  free  surface  reflections,  are  included.  The 


resulting  Green  functions  represent  the  complete  medium  response.  Since  the  wavenumber  response  is 
computed  as  a  function  of  frequency,  the  frequency-domain  Green  functions  used  in  equation  (2)  are 
obtained  directly. 

The  P  and  S-wave  velocity  models  estimated  for  The  Geysers  in  O’Connell  (1986)  from  a  P  and 
S-wave  progressive  velocity-hypocenter  inversion  are  used  in  Green  function  calculations  and  are  shown 
in  Figures  (1)  and  (2).  The  velocity  models  used  in  a  feasibility  study  of  moment  tensor  inversion  at 
The  Geysers  are  also  shown  in  Figures  (1)  and  (2).  The  initial  velocity  model  does  not  appear  to  be 
greatly  different  from  the  final  model.  However,  moment  tensor  inversions  with  the  initial  model  failed 
to  produce  the  correct  amplitude  pattern  of  P  and  S-wave  phases  on  any  components  of  ground  motion, 
in  contrast  to  the  results  obtained  using  the  final  velocity  models  (Figure  (7)).  The  progressive  P  and  S- 
wave  velocity-hypocenter  inversion  constituted  an  essential  component  of  the  moment  tensor  inversion 
process. 

Anelastic  attenuation  was  included  by  specifying  a  Q  model  for  The  Geysers  consistent  with  the 
results  of  Majer  and  McEvilly  (1979).  Values  of  50-100  were  used  for  Q P,  and  values  of  40-80  were 
used  for  05.  The  low  Q  values  were  used  near  the  free  surface  and  the  higher  values  used  in  the  pro¬ 
duction  zone.  These  low  Q  values  reflect  the  combined  effects  of  intrinsic  attenuation  and  scattering 
losses  and  correspond  to  effective  Q  values. 

Proper  phase  matching  of  observed  S-P  times  with  Green  function  S-P  times  is  important  to 
ensure  the  success  of  the  moment  tensor  inversions.  Earthquake  locations  estimated  in  O’Connell 
(1986)  were  used  to  define  initial  hypocenter-receiver  azimuths  and  distances.  Since  recording  station 
were  located  at  different  elevations,  the  predicted  S-P  times  for  initial  hypocenter-receiver  distance  did 
not  always  match  observed  S-P  times.  Hypocenter-receiver  distances  were  modified  so  as  to  produce 
correct  Green  function  S-P  times.  Hypocenter-  receiver  azimuths  are  preserved  but  takeoff  angles  are 
somewhat  different.  For  one  station  take-off-angles  were  altered  by  as  much  as  29°  for  one  event,  but 
for  most  stations  take-off-angles  were  not  changed  by  more  than  5°  -  15°.  Table  (1)  contains  the 
minimum  and  maximum  take-off-angle  error  ranges,  and  mean  take-off-angle  errors  and  their  standard 


deviations  for  all  events. 


VELOCITY  (KM/SEC) 


DEPTH  (KM) 


Ligure  1.  P  and  S-wave  velocity  models  used  for  a  feasibility  study  of  moment  tensor 
inversions  at  The  Geysers  (denoted  as  initial)  and  the  models  used  for  the  final  moment 
tensor  inversion  (denoted  as  final).  Note  that  the  initial  P  and  S-wave  velocity  models 
have  overall  velocity  gradients  similar  to  the  final  models. 
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DEPTH  (KM) 


Figure  2.  VyV,  initial  and  final  models  used  in  real  and  synthetic  waveform  moment 
tensor  inversions.  Note  that  the  assumption  of  constant  VyV,  used  in  the  initial  model  is 
clearly  in  error. 


Since  a  frequency  domain  inversion  is  used,  it  would  be  difficult  to  use  windows  about  certain 
phases  in  the  inversion.  Small  time  windows  about  the  first  P  and  S-wave  pulses  would  not  provide  the 
frequency  bandwidth  or  resolution  needed  to  reliably  estimate  the  moment  tensor.  Truncation  effects 
due  to  windowing  are  accentuated  for  short  time  windows.  Consequently,  complete  seismograms  were 
used  for  all  components  in  inverting  for  the  moment  tensor.  Ten  seconds  of  data  were  used  in  the 
inversions.  For  stations  that  had  shorter  records,  zeros  were  added  to  give  total  lengths  of  10  seconds. 

While  record  lengths  of  10  seconds  were  used  in  the  moment  tensor  inversions,  moment  tensor 
results  are  plotted  for  times  less  than  one  second.  This  was  done  because  source  durations  are  short, 
approximately  0.1  sec  and  because  the  Green  functions  do  not  contain  coda  wave  durations  as  long  as 
seen  in  the  observed  data.  Also,  some  components  of  the  observed  data  had  small  noise  glitches 
approximately  two  seconds  after  the  first  S-wave  arrival  and  these  glitches  contaminate  the  moment  ten¬ 
sor  time  functions  after  several  seconds.  Windows  of  less  that  one  second  were  used  to  detrend  the 
moment  rate  tensor.  Since  the  source  durations  of  these  microearthquakes  were  short,  this  approach  is 
reasonable. 


Take-off-angle  Errors 

Event 

Range 

Mean 

SD 

min 

max 

1 

-16° 

27° 

-1.2° 

11.5° 

2 

-16° 

-1° 

-8.1° 

5.5° 

-29° 

-15.6° 

9.8° 

Table  1.  Combined  P  and  S-wave  first  arrival  take-ofT-angle  errors  for  all  events,  defined  as  the 
difference  between  progressive  inversion  estimates  and  Green  function  estimates. 
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4.  Inversion  Results 


Moment  tensor  inversions  were  done  for  three  earthquakes  at  The  Geysers.  Two  of  the  cents 
were  shallow,  approximately  2  km  below  the  mean  station  elevation,  and  the  third  event  was  deeper, 
approximately  4  km  below  the  mean  station  elevation.  The  shallow  events  correspond  to  the  depth 
interval  where  the  strike-slip,  normal-slip,  reverse-slip  ambiguity  is  most  pronounced.  Gppenheimer 
(1986)  found  that  most  events  in  this  depth  range  had  strike-slip  solutions  with  a  smaller  number  of 
events  having  reverse-slip  mechanisms.  The  deeper  event  corresponds  to  the  depth  interval  where  nor¬ 
mal  faulting  mechanisms  predominate  (Oppenheimer,  1986). 

Event  and  station  locations  are  shown  in  Figure  (3).  Events  1  and  2  have  a  range  in  epicentral 
distances  representing  first  arrival  take-off  angle  coverage  over  45°,  but  also  have  a  large  azimuthal  gap 
in  station  coverage  to  the  south.  Conversely,  event  3  has  excellent  azimuthal  coverage,  but  very  limited 
take-off  angle  coverage.  All  event  inversions  were  full  rank,  with  the  condition  numbers  for  the  event  3 
inversion  about  twice  those  of  the  event  1  and  2  inversions  over  the  frequency  band  2-20  Hz.  The 
larger  condition  number  for  event  3  reflects  a  near  linear  dependency  due  to  the  limited  take-off  angle 
coverage.  Consequently,  the  event  3  inversion  is  most  likely  to  be  adversely  affected  by  errors  in  the 
data  or  Green  functions. 

Results  of  the  moment  tensor  inversion  are  displayed  in  the  following  manner.  The  orientations  of 
the  eigenvectors  of  DtJ  are  plotted  on  stereographic  lower  hemisphere  projections  along  with  available 
P-  wave  first  motion  data.  The  P-wave  first  motion  data  come  from  the  temporary  network  and  USGS 
stations.  In  order  to  obtain  as  many  first  motions  as  possible,  USGS  stations  outside  the  primary  pro¬ 
duction  zone  at  The  Geysers  were  used.  The  P-wave  velocity  model  estimated  in  O’Connell  (1986),  is 
not  completely  adequate  to  determine  azimuth  and  takeoff  angles  for  stations  outside  The  Geysers  for 
two  reasons.  Firstly,  P-wave  velocity  estimates  from  the  progressive  velocity-  hypocenter  inversion  are 
only  available  to  a  depth  of  4.0  km  and  more  distant  station  arrivals  correspond  to  rays  bottoming 
below  this  model  depth,  in  a  part  of  the  model  that  was  added  as  a  rough  estimate  of  the  velocity  struc¬ 
ture  below  4.0  km.  Secondly,  Eberhart-Phillips  (1986)  has  found  significant  lateral  variations  of  P-wave 
velocity  structure  outside  The  Geysers  production  zone,  so  azimuthal  estimates  may  be  in  error  for 
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Figure  3.  Map  showing  relation  of  earthquake  epicenters  (circles  with  event 
numbers  on  the  right)  to  stations  (A  and  V)  used  in  moment  tensor  inversion.  The  station 
denoted  by  the  V  was  not  used  for  Event  2.  Hypocentrai  depths  and  standard  errors  for 
events  1,  2,  and  3,  where  1.4  ±0.1,  1.6±0.2,  and  3.5±0.1  km,  respectively.  Ellipses 
represent  two-standard-deviation  epicentral  error  estimates.  The  solid  trace  is  Big  Sulfur 
Creek.  The  tine-dashed  lines  are  34.5  bar  contours  (Lipman  et  al.,  1978)  enclosing  areas 
of  pressure  decline  for  Lhe  year  1977  and  provide  a  rough  outline  of  The  Geysers’  pri¬ 
mary  steam  production  zones. 
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USGS  stations  well  outside  The  Geysers  due  to  lateral  refraction.  The  estimated  position  of  first 
motions  on  the  focal  sphere  of  the  distant  (>  20  km)  USGS  readings  may  have  substantial  uncertainties. 

Inversion  results  are  summarized  in  Table  (2)  and  Figures  (4-6)  for  events  1,  2,  and  3,  respec¬ 
tively.  The  estimated  principal  stress  orientations  for  these  three  Geysers  earthquakes  prove  to  be  quite 
stable  (Figures  (4-6)),  and  provide  results  consistent  with  observed  P-wave  first  motion  distributions. 
The  decompositions  of  the  estimated  moment  tensors  into  isotropic  and  deviatoric  components  show 
small  to  moderate  departures  from  a  single  double-couple  source  for  events  1  and  2  and  large  ones  for 
event  3  (Table  (2)).  The  sometimes  large  isotropic  component  for  event  3  can  be  partially  explained  by 
instability  of  the  Ma  component  due  to  the  aforementioned  take-off  angle  problem.  However,  it  is  not 
immediately  clear  how  to  interpret  the  apparent  source  mechanism  complexity  in  light  of  potential 
sources  of  error  in  the  Green  functions,  such  as  the  take-off  angle  errors  listed  in  Table  (1),  source 
mislocation,  and  velocity  model  errors.  In  an  effort  to  quantify  the  effects  of  Green  function  a  series  of 
synthetic  tests  were  done  and  are  described  in  the  next  section.  Consequently,  we  defer  further  interpre¬ 
tation  of  the  moment  tensor  estimates  to  a  latter  section  in  order  to  incorporate  the  synthetic  test  results. 
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Figure  4.  Lower-hemisphere  equal-area  plot  of  P-wave  first  motions,  time  varying 
principal  moment  axes  (0.5  second  duration)  estimated  from  the  moment  tensor  inver¬ 
sion,  and  double-couple  fault  plane  solution  for  event  1.  Compressional  and  dilata- 
tional  P-wave  first  motions  are  plotted  as  (C  and  +)  and  (D  and  -),  respectively.  The 
tension  axis  is  denoted  by  the  large  (T)  and  its  time  history  is  shown  as  line  segments 
with  arrowheads  porn  Ling  toward  the  next  point  in  time.  The  compression  axis  is  denote 
by  a  large  (P)  and  the  intermediate  axis  is  denoted  by  a  large  (I).  The  solid  lines  are  no¬ 
dal  planes  drawn  to  satisfy  the  first  motion  data.  The  inconsistent  compression  in  the 
upper  right  quadrant  corresponds  to  a  distant  station  and  its  position  has  large  uncertain¬ 
ties  due  to  lateral  refraction  effects. 
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Moment  Tensor  Inversion  Summary 

Event 

Sta. 

Comp. 

Mc 

K 

Isotropic 

CLVD 

T 

o 

H 

min 

max 

min 

max 

(sec) 

1 

7 

7 

14 

1.6 

0.2 

0 

15 

10 

K9 

m 

2 

6 

6 

11 

1.8 

0.8 

10 

25 

10 

rj 

3 

7 

7 

LJO 

2.3 

3.0 

10 

55 

0 

Efl 

Table  2.  Moment  inversion  information  by  event  number  showing  the  number  of  stations  and  ground 
motion  components  used,  V  for  vertical  and  H  for  horizontal.  Mc  is  U.S.G.S  coda  magnitude,  Ma  is 
scalar  moment  in  units  of  1020  dyne-cm.  The  last  5  columns  represent  minimum  and  maximum  isotropic 
percentages  of  the  total  moment,  minimum  and  maximum  CLVD  percentages  of  the  deviaioric  moment, 
and  t  is  the  moment  time  duration  used  for  these  estimates. 
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5.  Synthetic  Tests 


The  formulation  for  solving  for  the  moment  tensor  presented  in  equations  (1-2)  assumed  that 
Green  function  errors  were  insignificant.  Using  matrix  notation,  a  solution  for  the  moment  tensor  ele¬ 
ments,  m.  is  written 

m  =  C'u  (5) 

where  G"1  is  the  inverse  Green  function  and  u  are  the  measured  ground  motion  data.  Moment  tensor 
inversion  approaches  which  assume  that  errors  only  occur  in  the  ground  motion  data,  u,  will  underesti¬ 
mate  errors  in  moment  tensor  estimates  and  their  eigenvalue-eigenvector  decompositions.  A  more  realis¬ 
tic  representation  is  to  rewrite  (5)  in  the  form 

M=/(G-‘,u)  (6) 

to  emphasize  that  G  is  also  an  unknown  to  some  degree.  We  investigate  the  mapping  of  uncertainties  in 
G  into  estimates  of  M  by  determining  estimates  of  M  for  a  suite  of  extremal  values  of  G  using  (5). 
The  effects  of  Green  function  errors  can  be  assessed  from  the  resulting  variation  of  M  estimates  thus 
obtained. 

For  microearthquakes  the  primary  sources  of  Green  function  error  are  earthquake  mislocation  and 
incorrect  velocity  structure.  Source  mislocations  effects  are  significant  because  stations  are  close  to  epi¬ 
centers,  producing  Green  function  errors  due  to  inconect  distances,  azimuths,  and  take-off  angles. 
Incorrect  velocity  models  produce  amplitude  and  take-off  angle  errors  in  calculated  Green  functions. 

The  synthetic  tests  focus  on  source  mislocation  effects,  although  a  preliminary  effort  to  appraise 
velocity  model  error  effects  is  included.  There  are  two  reasons  for  focusing  on  mislocation  effects. 
Firstly,  realistic  estimates  of  hypocentral  uncertainties  for  the  three  Geysers  events  are  available  from  a 
progressive  velocity-hypocenter  inversion  (Eberhart-Phillips,  1986;  O’Connell,  1986)  and  a  nonlinear 
hypocentral  error  appraisal  (O’Connell,  1986).  Secondly,  although  velocity  models  have  been  deter¬ 
mined  for  the  Geysers  (Eberhart-Phillips,  1986;  O’Connell,  1986),  the  intensive  computational  require¬ 
ments  of  completely  appraising  velocity  model  errors  is  beyond  the  scope  of  the  present  investigation. 
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The  95%  confidence  ellipses  from  O’Connell  (1986)  for  the  epicentral  locations  of  all  three 
events  are  shown  in  Figure  (3)  and  95%  confidence  hypocentral  depth  uncertainties  are  presented  in 
Table  (3).  Eberhart-Phillips  (1986)  noted  that  a  one-dimensional  velocity-hypocenter  inversion  with  sta¬ 
tion  corrections  (the  method  used  in  O’Connell,  1986)  produce  epicentral  estimates  close  to  those 
obtained  using  a  three-dimensional  velocity-hypocenter  inversion.  However,  a  possibility  of  a  sys¬ 
tematic  epicentral  bias  exists;  epicentral  uncertainties  may  be  a  factor  of  two  larger.  Consequently,  syn¬ 
thetic  error  estimates  reflect  the  minimum  error  that  could  be  attributed  to  hypocentral  mislocation  at  a 
95%  confidence  level. 

The  double-couple  mechanisms  determined  from  the  moment  tensor  inversions  and  first-motion 
data  (solid  nodal  lines  in  Figures  (4-6))  were  used  to  define  the  "true"  source  mechanisms  of  the  three 
events  in  the  synthetic  tests.  The  "true”  synthetic  waveform  data  were  generated  using  these  source 
mechanisms,  the  "final"  velocity  models  in  Figure  (1),  and  the  "correct"  hypocenters.  Simple  double- 
couple  source  mechanisms  were  used  to  determine  what  percentage  of  spurious  isotropic  and  CLVD 
components  could  be  produced  simply  by  using  incorrect  earthquake  locations  and/or  incorrect  velocity 
models.  The  same  stations  and  ground  motion  components  used  in  the  real  data  inversions  were  used 
in  the  synthetic  inversions. 

To  determine  the  effects  of  mislocation  errors,  moment  tensor  inversions  were  done  using  Green 
functions  calculated  using  hypocentral  positions  on  the  95%  confidence  error  ellipsoid.  The  combined 
effects  of  mislocation  and  velocity  model  errors  were  investigated  by  using  the  initial  velocity  models 
in  Figure  (1)  to  calculate  the  mislocation  Green  functions.  This  approach  was  implemented  using  the 
WKBJ  method  (Dey-Sarkar  and  Chapman,  1978)  to  calculate  the  Green  functions  since  all  potential  first 
arrivals  represented  upgoing  ray  paths.  Because  many  values  of  G  must  be  computed,  the  WKBJ 
method  was  chosen  for  its  combination  of  accuracy  and  computational  efficiency.  To  further  reduce  the 
computational  burden,  known  source  time  functions  were  assumed,  allowing  time  domain  inversions 
for  the  static  moment  tensor.  Eight  point  lime  windows  beginning  at  the  first  P  and  S-wave  arrivals 
were  used  in  the  inversions.  In  contrast  to  inversions  with  real  waveform  data,  the  waveform  data  are 
error  free  in  the  synthetic  tests. 
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Specific  results  of  particular  synthetic  tests  cannot  be  generalized  because  they  depend  on  source- 
receiver  geometry  and  source  mechanism.  At  a  fixed  source -receiver  geometry  changing  the  source 
mechanism  will  produce  a  different  radiation  pattern  relative  to  the  station  distribution  on  the  focal 
sphere.  Similarly,  for  a  fixed  source  mechanism  one  station  distribution  can  be  more  prone  to  errors 
than  another.  Consequently,  moment  tensor  error  characteristics  need  to  be  investigated  on  a  event  by 
event  basis. 


Synthetic  Results 


GF  err. 

Isotropic 

CLVD 

min 

max 

min 

max 

misloc. 

0.1 

0.7 

0.0 

+vel. 

0.2 

9.5 

10.6 

misloc. 

0.0 

■a 

0.5 

7.7 

+vel. 

9.5 

H 

18.6 

26.5 

7.7 

3.9 


Table  3.  Synthetic  moment  tensor  inversion  test  results  by  event  number  showing  ranges  of  spurious 
isotropic  and  CLVD  component  percentages  and  maximum  P  and  T  axis  errors.  For  each  event,  results 
arc  listed  for  the  cases  of  mislocation  and  combined  mislocation  and  velocity  model  Green  function 
errors.  Maximum  depth  mislocations  are  listed  as  A?  (see  Figure  (3)  for  epicentral  mislocations). 


Results  of  the  synthetic  tests  are  summarized  in  Table  (3).  Essentially  the  same  results  were 
obtained  when  a  zero  isotropic  constraint  was  used  and  are  therefore  omitted.  The  results  listed  in 
Table  (3)  represent  a  lower  error  limit  for  the  source-receiver  geometry  considered.  Overall  errors  for 
real  data  moment  tensors  inversions  would  be  larger  due  to  the  combined  effects  of  Green  function 
errors  and  statistical  errors  in  the  data. 

Errors  in  P  and  T  axis  positions  are  relatively  smaller  than  errors  in  the  physical  source  decompo¬ 
sition  into  isotropic,  CLVD,  and  simple  double  couple.  The  results  of  these  and  other  synthetic  tests 
indicate  that  principal  moment  axis  orientations  are  much  less  sensitive  to  mislocation  and  velocity 
model  Green  function  errors  than  the  physical  source  decomposition  into  isotropic,  CLVD,  and  simple 
double-couple  components.  This  result  make  intuitive  sense  because  a  variety  of  physical  source  com¬ 
ponents  can  be  assembled  that  are  consistent  with  a  particular  set  of  principal  moment  orientations  and 
produce  certain  level  of  misfit  to  the  observed  data.  For  instance,  adding  small  spurious  isotropic 
source  components  might  compensate  for  systematic  underestimation  of  P-wave  amplitudes  due  to  a 
systematic  Green  function  error.  Yet,  the  principal  moment  axis  orientations  will  be  unaffected. 

The  waveform  misfits  produced  by  using  incorrect  velocity  models  and  hypoccntcrs  to  calculate 
Green  functions  were  small.  The  erroneous  isotropic  and  CLVD  components  effectively  compensated 
for  Green  function  errors  to  yield  small  waveform  misfits.  For  the  source-receiver  geometries  of  these 
events,  small  waveform  misfits  alone  were  not  a  reliable  indicator  that  moment  tensor  errors  were  small. 

6.  Discussion 

From  Table  (3)  it  is  clear  that  moment  tensor  decomposition  errors  due  to  mislocation  alone  are 
small  for  these  three  events,  although  small  spurious  CLVD  percentages  were  obtained  for  event  2  and 
3.  When  velocity  model  errors  are  combined  with  locations  errors,  the  spurious  isotropic  and  CLVD 
components  are  large  enough  to  explain  the  sizes  of  the  isotropic  and  CLVD  components  obtained  for 
events  1  and  2  using  real  data.  A  portion  of  the  CLVD  component  for  all  events  may  be  caused  by  dis¬ 
tortion  of  a  dominantly  quadrapole  radiation  pattern  due  to  lateral  velocity  heterogeneity. 


The  larger  isotropic  and  CLVD  components  estimated  for  event  3  may  reflect  poorer  velocity 
resolution  in  the  progressive-velocity  hypocenter  inversion  since  the  source  depth  is  located  in  the 
deepest,  least  resolvable  portion  of  the  estimated  velocity  model.  Alternatively,  there  is  no  reason  to 
assume  a  priori  that  the  larger  isotropic  and  CLVD  component  estimate  for  event  3  might  not  reflect  a 
true  source  property.  The  deeper  seismicity  at  The  Geysers  appears  to  delimit  the  lower  boundaries  of 
the  steam  field  and  appears  to  have  predominantly  normal-faulting  mechanisms  (Oppenheimer,  1986). 
The  large  CLVD  component  for  event  may  indicate  tensile  failure  is  occurring  at  the  base  of  the  steam 
field.  Tensile  failure  has  been  inferred  at  another  geothermal  field  by  Foulger  and  Long  (1984).  How¬ 
ever,  Foulger  and  Long  clearly  delineated  a  strong  non-quadrapole  radiation  pattern  from  events  with  a 
comparable  distribution  of  first-motion  data,  but  the  first-motion  data  of  event  3  is  compatible  with  a 
quadrapole  radiation  pattern  (Figure  (6)).  Thus,  known  and  potential  errors  in  the  assumed  Green  func¬ 
tions,  and  the  unfavorable  source-receiver  geometry  for  event  3  probably  produced  substantial  spurious 
isotropic  and  CLVD  components.  Further  refinements  of  the  moment  tensor  inversion  approach  and 
improved  Green  function  calculations  will  be  required  to  determine  the  significance  of  the  isotropic  and 
CLVD  components  of  event  3. 

Since  none  of  the  deviatoric  moment  tensors  corresponds  to  a  single  double-couple  for  the  entire 
moment  duration,  the  P-wave  nodal  surfaces  will  not  correspond  to  orthogonal  planes  for  the  entire 
moment  duration.  However,  results  of  the  synthetic  tests  indicate  that  significant  deviations  from  single 
double-couple  source  mechanisms  cannot  be  resolved  for  these  three  events.  Further,  there  is  no  evi¬ 
dence  from  the  P-wave  first  motion  distributions  for  significant  deviations  from  a  quadrapole  radiation 
pattern.  Consequently,  in  the  following  analysis  we  assume  that  a  simple  double-couple  source  mechan¬ 
ism  is  an  appropriate  approximation  to  interpret  the  deviatoric  moment  tensor. 

The  P-wave  first  motions  for  events  1  and  3  provide  sufficient  constraints  to  unambiguously 
resolve  the  focal  mechanism  (Figures  (4)  and  (6)).  The  estimated  principal  moments  orientations  for 
events  1  and  3  arc  in  excellent  agreement  with  the  first  motion  focal  mechanism  solutions.  These 
results  bear  out  the  synthetic  test  conclusions  that  principal  moment  orientation  estimates  are  robust 
with  respect  to  Green  function  errors.  Note  that,  if  the  north-most  and  west-most  dilatations  were 
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unavailable  in  Figure  (4),  a  normal  faulting  mechanism  would  be  compatible  with  the  remaining  P-wave 
first  motions.  The  moment  tensor  solution  provides  sufficient  information  to  resolve  the  correct  focal 
mechanism  and  required  waveform  information  from  only  7  stations.  In  contrast,  21  P-wave  first 
motions  were  required  to  constrain  the  nodal  planes  in  Figure  (4). 

For  event  2  the  principal  moment  positions  are  used  as  constraints  to  construct  the  "best"  focal 
mechanism,  using  the  P-wave  first  motions  as  the  additional  constraints  (Figure  (5).  The  fault  plane 
solutions  shown  in  Figure  (5)  demonstrate  that  a  wide  range  of  focal  mechanisms  are  consistent  with 
the  P-wave  first  motion  data.  Solutions  ranging  from  nearly  pure  normal-slip  to  pure  strike-slip  are 
consistent  with  the  first  motion  data  (Table(4)).  The  moment  tensor  solution  constrains  the  solution  to 
be  dominantly  strike-slip.  The  orientations  of  the  principal  moment  axes  are  stable  even  though  there  is 
a  large  azimuthal  gap  in  station  coverage. 

Event  2  is  the  type  of  earthquake  that  made  Oppenheimer’s  (1986)  reduction  of  his  fault  plane 
solution  data  difficult.  The  P-wave  first  motions  are  consistent  with  both  strike-slip  mechanisms  (postu¬ 
lated  shallow  event  mechanism)  and  normal-slip  mechanisms  (postulated  deep  event  mechanism).  He 
would  have  been  forced  to  discard  this  event  because  constraints  are  absent  with  respect  to  the  depth- 
dependent-focal-mechanism  hypothesis.  The  moment  tensor  waveform  inversion  method  provided 
sufficient  resolution  to  distinguish  the  appropriate  focal  mechanism  using  data  from  a  subset  of  stations. 

The  fault  plane  solution  in  Figure  (6)  is  constrained  by  three  first  motions,  the  two  west-most 
dilatations,  and  the  northeast  compression  shown  as  (+).  If  these  3  first  motions  were  unavailable,  an 
almost  pure  strike-slip  mechanisms  would  fit  the  first  motion  data.  These  three  first  motions  would  not 
be  available  for  a  USGS  solution,  since  the  dilatations  represent  temporary  station  readings,  and  the  (+) 
is  an  ambiguous  reading  from  a  distant  station.  The  full  first-motion  data  set  confirms  the  moment  ten¬ 
sor  P  and  T  axis  estimates  for  event  3  and  demonstrates  that  comparable  focal  mechanism  resolution 
can  be  achieved  with  the  waveform  inversion  using  data  from  7  stations,  instead  of  the  22  first-motions 
required  to  constrain  the  focal  mechanism  in  Figure  (6). 
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P  and  T  Axis  Orientations 


First  Motion 

Waveform 

P 

T 

_ l _ 

T 

Event 

trend 

plunge 

trend 

plunge 

trend 

plunge 

2131638 

47° 

2° 

133° 

14° 

50°, 62° 

145°,150° 

2181937 

59° 

0° 

149° 

5° 

53°,64° 

2°  2° 

142°, 154° 

-4°, 11° 

18° 

72° 

132° 

11° 

2200908 

355° 

52° 

92° 

5° 

338°,358° 

-1°,14° 

Table  4.  Comparison  of  P  and  T  axis  orientations  estimated  from  first-motion  fault  plane  solutions  and 
waveform  inversions.  The  range  of  time  variations  of  orientations  of  the  principal  moment  axes  are 
listed  for  the  waveform  moment  tensor  estimates.  Event  2  has  an  ambiguous  first-motion  fault  plane 
solution  so  two  possible  first  motion  solutions  are  listed. 
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7.  Comparison  of  Observed  and  Predicted  Seismograms 


The  earliest  waveform  moment  tensor  inversion  attempts  used  the  constant  velocity  layer  parame¬ 
terization  model  for  The  Geysers  of  Eberhart-Phillips  and  Oppenheimer  (1984).  Waveform  amplitude 
predictions  for  P  and  S-wave  phases  were  very  poor  even  when  the  best  moment  tensors  estimates 
obtained  with  the  final  model  were  used  The  velocity  discontinuities  of  the  relatively  thick  (1.0  km) 
constant  velocity  layers  produced  phases  arrivals  that  were  more  representative  of  parameterization  than 
actual  earth  structure.  The  initial  model  in  Figure  (1)  was  produced  as  a  linear  velocity  gradient 
modification  and  improved  waveform  amplitude  predictions,  but  misfits  were  clearly  unacceptable. 
However,  when  the  final  models  (Figure  (1-2))  of  a  progressive  vclocity-hypocenter  inversion 
(O'Connell,  1986)  were  used,  excellent  agreement  was  obtained  between  relative  P  and  S-wave  first- 
arrival  amplitudes  for  all  inversions  (Figure  (7)).  The  primary  differences  between  the  initial  and  final 
models  are  reflected  in  Vp/V,  (Figure  (2))  ar.J  larger  velocity  gradients  near  a  depth  of  one  kilometer 
for  the  final  model.  .‘gTeemeni  between  observed  and  predicted  seismograms  using  the  final  velocity 
model  is  very  good  for  both  P  and  S-wave  first  arrivals. 

The  initial  models  shown  in  Figure  (1)  arc  incorrect  for  The  Geysers.  They  would  be  much 
closer  to  the  estimated  final  models  than  a  small  set  of  constant  velocity  layers,  a  parameterization  com¬ 
monly  employed  to  calculate  Green  functions  for  waveform  studies  of  local  events.  The  constant  V,/V5 
assumption  made  for  the  initial  model  in  Figure  (2)  was  clearly  wrong  for  The  Geysers,  but  also  reflects 
an  assumption  often  made  in  waveform  studies. 

Saikia  and  Herrmann  (1980)  suggested  that  unsatisfactory  S-wave  synthetic  fils  from  moment  ten¬ 
sor  waveform  inversions  for  two  Arkansas  microearthquakc  were  probably  a  reflection  of  inaccurate  S- 
wavc  velocity  models.  Our  experiences  arc  consistent  with  their  conclusions  because  waveform  model¬ 
ing  deficiencies  resulting  from  using  the  constant  velocity  layer  and  imual  model  parameter) zations 
were  very  similar  to  those  of  Saikia  and  Herrmann  (1986).  These  results  reveal  the  importance  of  using 
realistic  velocity  models  for  microearthquakc  moment  tensor  inversions. 
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!y.  Note  that  the  relative  amplitudes  of  i he  direct  P  and  S-wavs  arrivals  on  com¬ 
ponents  of  ground  mouon  are  faithfully  reproduced  in  the  synthetic  seismograms. 
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8.  Summary  and  Conclusions 


Frequency  domain  moment  tensor  inversions  were  done  for  three  microcarthquakes  from  The 
Geysers  geothermal  field.  The  available  first-motion  data  placed  strong  constraints  on  two  microearth¬ 
quake  focal  mechanisms.  Waveform  principal  moment  orientation  estimates  agreed  with  P-wave  first 
motions  focal  mechanisms  for  these  two  events  and  constrained  the  focal  mechanism  of  the  event  with 
ambiguous  first-motion  solutions.  The  moment  tensor  estimates  of  principal  moment  orientations  were 
obtained  using  far  fewer  stations  than  required  for  first-motion  focal  mechanisms  solutions. 

The  three  focal  mechanisms  obtained  here  support  Oppenheimer's  (1986)  hypothesis  that  focal 
mechanisms  arc  a  function  of  depth  at  The  Geysers.  Specifically,  strike-slip  focal  mechanisms  were 
obtained  for  the  two  shallow  events  and  a  predominantly  normal-slip  focal  mechanism  was  obtained  for 
the  deep  event,  as  his  model  predicts.  The  orientation  of  the  T  for  the  shallow  events  (events  1  and  2 
in  Table  (1))  is  rotated  40’  to  50°  clockwise,  with  respect  to  his  estimate  of  105°  as  the  azimuth  of  least 
compressive  stress  for  The  Geysers.  Since  only  two  events  are  available  for  comparison  the  differences 
are  probably  not  significant  The  potential  of  moment  tensor  inversions  to  provide  well  constrained 
principal  moment  orientations  for  individual  events  may  make  improve  our  understanding  of  the  rela¬ 
tionship  between  seismicity,  steam  production,  and  tectonic  processes  at  The  Geysers. 

The  results  obtained  here  were  dependent  on  reliable  estimates  of  velocity  structure  so  as  to 
minimize  errors  in  calculated  Green  functions.  The  velocity  model  determined  in  O’Connell  (1986), 
when  used  in  the  moment  tensor  inversions,  produced  the  correct  ratio  of  P  and  S-wave  amplitudes  on 
all  components  of  ground  motion  (Figure  (7)).  This  result  strengthens  the  arguments  in  O’Connell 
(1986)  that  the  estimated  models  are  good  one -dimensional  representations  of  the  velocity  structure  at 
The  Geysers. 

The  satisfactory  results  obtained  here  are  a  direct  consequence  of  using  P  and  S-wave  data 
together  to  estimate  velocity  structure,  hypocentcr  locations,  and  moment  tensors.  Although  moment 
tensor  inversions  were  not  done  to  compare  the  effects  of  using  just  P-wavc  data  to  those  using  both  P 
and  S-wave  data.  Stump  and  Johnson  (1977)  found  that  inversions  that  just  used  P-wave  maximum 
amplitudes  were  not  as  well  conditioned  as  inversions  using  complete  seismograms.  It  was  clear  from 
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the  synthetic  fits  to  the  data,  that  the  S-wave  phases  significantly  constrained  the  moment  tensor  esti¬ 


mates. 

Further  improvements  in  assessing  the  effects  of  Green  function  errors  on  moment  tensor  esti¬ 
mates  are  needed.  The  synthetic  tests  presented  here  were  used  to  roughly  quantify  Green  function 
error  effects.  An  obvious  next  step  is  to  incorporate  the  methods  used  in  the  synthetic  tests  directly  into 
real  data  moment  tensor  inversions.  Modifications  of  the  present  moment  tensor  inversion  procedure 
will  allow  elimination  of  take-off  angle  errors  due  to  large  station  elevation  differences.  These 
improvements  are  required  to  resolve  the  significance  of  moderate  amounts  of  isotropic  and  CLVD 
components  contained  in  the  microearthquake  moment  tensor  estimates. 

Conditions  at  The  Geysers  are  similar  to  many  geologically  important  areas  of  microseismicity. 
Consequently,  The  Geysers  represents  a  good  test  of  the  feasibility  of  doing  frequency-domain  moment 
tensor  inversions  in  geologically  complex  areas.  Topographic  variations  are  large,  reflected  in  the  0.6 
km  variation  of  station  elevations,  near  surface  velocity  variations  arc  profound,  and  significant  lateral 
velocity  heterogeneities  exist.  These  factors  were  not  included  in  the  moment  tensor  inversions,  and 
useful  results  were  obtained  despite  rather  unfavorable  station  recording  geometries  for  all  three  events 
considered.  The  real  and  synthetic  data  inversions  demonstrate  that,  with  the  proper  attention  to  deter¬ 
mining  realistic  velocity  structure  and  event  locations,  moment  tensor  estimates  for  microearthquake 
sources  can  be  robust  and  provide  a  means  to  resolve  source  mechanisms  using  data  from  a  relatively 
small  number  of  stations. 
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